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During  the  first  year  of  this  effort,  two  relatedissues  were  investigated:  (1)  establishing  the 
validity  of  the  thin  flame  model  when  used  to  compote  flow-combustion  interactions  in  a 
turbulent  shear  layer;  (2)  developing  an  efficient  methodology  to  compute  the  unsteady  strained 
flame  structure  when  the  flame  thickness  is  much  smaller  than  the  flow  scale.  In  the  first  effort, 
the  transport  element  method  was  applied  to  compute  (a)  a  reacting  flow  in  which  combustion 
proceeds  according  to  a  single-step,  temperature  dependent  Arrhenius  reaction,  and  (b)  a  mixing- 
limited  model  in  which  Schvab-Zeldovich  variables  are  used  to  obtain  the  infinite  speed 
chemistry  results.  The  results  of  both  computations  showed  that,  at  high  Damkohler  numbers, 
while  there  is  a  small  error  in  the  prediction  of  the  total  burning  rate  using  the  second  approach, 
the  second  model  accurately  estimates  the  effect  of  combustion  on  the  flow  dynamics  in  terms  of 
volumetric  expansion  and  vorticity  generation.  Work  on  the  second  project  resulted  in  a  more 
efficient  model  to  compute  the  flame  structure  under  conditions  of  unsteady  strain.  The 
computational  model  is  based  on  a  series  of  mathematical  transformations  which  reduce  the 
governing  equations  to  time-dependent  reaction-diffusion  equations. 
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During  the  first  year  of  this  effort,  we  focused  on  two  related  problems:  (1)  establishing 
the  validity  of  the  thin  flame  model  when  used  to  compute  flow-combustion  interactions  in  a 
turbulent  shear  layer;  (2)  developing  an  efficient  methodology  to  compute  the  unsteady  strained 
flame  structure  when  Ae  flame  diickness  is  much  smaller  than  the  flow  scale.  In  the  first  effort, 
the  transport  element  method  was  applied  to  compute  (a)  a  reacting  flow  in  which  combustion 
proceeds  according  to  a  single-step,  temperature  dependent  Arrhenius  reaction,  and  (b)  a  mixing- 
limited  model  in  which  Schvab-Zeldovich  variables  are  used  to  obtain  the  infinite  spe^ 
chemistry  results.  The  results  of  both  computations  showed  that,  at  high  Damkohler  numbers, 
while  there  is  a  small  error  in  the  prediction  of  the  total  burning  rate  using  the  second  approach, 
the  second  model  accurately  estimates  the  effect  of  combustion  on  the  flow  dynamics  in  terms  of 
volumetric  expansion  and  vorticity  generation.  Thus,  implementing  a  detailed  flame  sheet  model 
using  a  flowfield  decomposition  should  yield  an  accurate  simulation  for  the  flowfield  while 
maintaining  the  overall  computational  requirements  below  what  is  needed  in  the  original 
analysis.  Work  on  the  second  project  has  resulted  in  a  more  efficient  model  to  compute  the 
flame  structure  under  conditions  of  unsteady  strain.  The  computational  model  is  based  on  a 
series  of  mathematical  transformations  which  reduce  the  governing  equa  ons  to  time-dependent 
reaction-diffusion  equations.  Model  results  have  been  u^  to  determine  the  flame  characteristic 
response  time  and  the  effect  of  strain-Lewis  number  on  the  burning  velocity.  The  model  is 
currently  being  extended  to  multistep  chemistry. 
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I.  OBJECTIVES 

The  objectives  of  this  work  are: 

(1)  To  develop,  using  methods  of  asymptotic  expansion,  a  computational  framework  for 
the  simulation  of  turbulent  combustion  phenomena  by  deriving  equations  which  govern  the  flow 
and  combustion  in  different,  physically-distinct  regions  in  the  domain.  This  will  be  independent 
of  the  numerical  method,  or  methods,  used  to  integrate  these  equations  and  hence  will  have 
applications  beyond  vortex  simulations; 

(2)  To  develop  subscale  fundamentally  based  models  which  can  be  used  to  obtain  large 
eddy  simulations  using  vortex  methods.  Subscale  models,  in  this  case,  will  be  obtained  through 
the  application  of  the  renormalization  group  theory  to  the  equations  governing  vorticity 
stretching  and  tilting  on  scales  lower  than  numerically  resolvable  scales.  The  application  of 
RNG  in  physical  vorticity  space  is  motivated  by  our  solutions  which  show  that  the  properties  of 
vortex  lines  at  the  small  scales  follow  closely  the  RNG  predictions. 

(3)  To  develop  fundamentally  based  flame  structure  models  using  the  equations  obtained 
in  (1)  where  flow  combustion  interactions  are  represented  by  the  effect  of  the  time-dependent 
stretch  exerted  by  the  outer  flow  on  the  flame  structure.  One  important  respect  of  these  models 
will  be  the  incorporation  of  multistep  chemical  kinetics  algorithms  to  accurately  capture  flow 
combustion  interactions  when  it  is  dominated  by  radical  concentration  and  diffusion; 

(4)  To  modify  the  transport  element  method;  a  Lagrangian  scheme  which  we  developed 
to  simulate  reacting  species  transport,  to  act  as  a  “coupling”  algorithm  between  the  solutions 
obtained  for  the  outer  flow  and  the  inner  flow. 

The  developed  methodology  will  be  applied  to  study  mixing  and  combustion  in  reacting 
shear  flow  at  high  Reynolds  numbers.  Effect  of  strong  density  variations,  high  heat  release  rates 
at  elevated  Damkohler  numbers,  and  high  Mach  numbers  will  also  be  investigated. 
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II.  PERSONNEL 

During  the  period  of  1992-1993,  funding  was  used  to  support  the  work  of  the  following 
students: 

(1)  Marios  Soteriou  who  completed  his  Ph.D.  thesis  in  June  1993  and  stayed  on  as  a 
postdoctor. 

(2)  Van  Luu,  a  Ph.D.  student  working  on  distributed  reaction  zone  combustion  model. 

(3)  Constantin  Petrov  who  completed  his  master’s  work  in  June  1993  and  stayed  on  as  a 
Ph.D.  student  working  on  the  thin  flame  combustion  model. 

in.  WQRK.STATII.S 

m.l  Flow  Combustion  Interaction  Modeling  Using  a  Flame  Sheet  Representation. 

A  numerical  methodology  has  been  introduced  to  enable  the  study  of  a  post-transitional 
spacially  developing  exothermically  reacting  shear  layer  over  a  substantial  range  of  the 
governing  parameters.  The  Transport  Element  Method,  commonly  used  in  the  simulation  of 
non-reacting  flows  is  further  developed  to  accommodate  an  exothermically  reacting  flowfield. 
The  scheme  is  Lagrangian,  grid-free  and  adaptive  and  solves  the  unaveraged,  time  dependent  and 
coupled  scalar  transport-reaction  and  Navier-Stokes  equations  respectively,  in  their  scalar- 
gradient  and  vorticity  forms.  It  exploits  the  Shvab-Zeldovich  formulation  to  provide  solutions 
for  both  moderately  fast  and  infinitely  fast  reactions.  For  finite  reaction  speeds,  Arrhenius 
kinetics  are  used. 

Numerical  results  were  used  in  a  preliminary  study  of  the  effects  combustion 
exothermicity  on  both  the  flow  and  scalar  fields.  We  found  that  the  externally  forced  flowfield  is 
substantially  modified  in  the  presence  of  combustion  and  a  reduced  growth  is  experienced 
mainly  due  to  an  alteration  of  the  mechanism  by  which  the  vortical  structures,  which  characterize 
the  flow,  interact.  This,  together  with  the  decreased  density  within  the  mixing  region  leads  to 
decreased  efficiency  of  combustion. 
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The  simpler,  computationally  less  demanding  inflnite  reaction  speed  model  was  found  to 
be  an  effective  model  of  the  detailed  finite  reaction  problem  as  long  as  the  reactions  involved  in 
the  latter  are  fast  compared  to  the  flow,  the  reaction  zone  is  thin,  and  abrupt  transient  phenomena 
such  as  quenching  are  avoided.  Detail  of  this  study  are  presented  in  Appendix  I  of  this  report. 


The  transient  response  of  a  flame  to  an  unsteady  strain  rate  was  analyzed  using  a  series  of 
mathematical  transformations.  Initially,  the  flame  was^located  in  the  stagnation  plane.  At  time 
t=0  the  unsteady  strain  rate  was  applied.  Two  basic  patterns  of  the  strain  rate  are  considered:  a 
step  function,  and  a  sinusoidal  wave.  /The  flame  response  was  characterized  by  two  parameters: 
(1)  the  burning  velocity  and,  (2)  the  flame  location.  The  influence  of  the  Lewis  number  and  the 
strain  rate  on  these  two  parameters  and  on  the  relaxation  time  was  investigated.  For  unity 
Lewis  number,  and  within  the  range  of  compressive  strains  characteristic  for  a  turbulent  jet 
flow  (200  -  500  1/sec),  the  shapes  of  the  temperature  and  mass  fraction  profiles  remain  almost 
unchanged.  This  leads  to  a  burning  velocity  which  is  only  slightly  dependent  on  the  strain  rate. 

In  the  case  of  non-unity  Lewis  number  the  profiles  are  significantly  altered  even  by  relatively 
weak  strain  rates.  This  happens  due  to  the  interaction  of  unbalanced  heat  and  concentration 
diffusion  and  convection.  The  changes  in  the  profiles  produce  significant  variation  in  the 
burning  velocity.  Some  analytical  results  are  obtained  for  the  relaxation  time  of  flame  as  a 
function  of  the  strain  rate  and  thermo-chemical  parameters.  The  analytical  derivations  are  based 
on  the  application  of  the  integral  approach  in  the  transformed  domain.  In  order  to  investigate  the 
receptivity  characteristics  of  flame,  the  periodic  strain  rate  is  applied.  Over  a  wide  range  of 
frequencies,  flame  demonstrates  periodic  response.  The  average  value  of  the  burning  velocity  is 
very  close  to  the  burning  velocity  of  a  flame  under  the  average  strain;  phase  shift  between  the 
burning  velocity  and  strain  rate  fluctuations  is  approximately  constant  and  equal  to  -1 .3;r  .  This 
pattern  is  violated  only  when  the  period  of  the  strain  rate  oscillations  is  much  lower  than  the 
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diffusion  time  scale.  The  analysis  of  the  flame  response  suggests  that  the  steady  state 
assumption  can  be  used  with  a  reasonable  accuracy  in  a  flamelet  modeling,  although  the  phase 
shift  should  be  taken  into  account.  Extinction  strain  rate  which  we  define  as  the  strain  rate  when 
the  steady  state  flame  location  crosses  the  stagnation  plane  for  the  first  time,  is  an  exponential 
function  of  the  heat  release.  This  suggests  that,  in  order  to  get  adequate  values  of  the  extinction 
strain  using  a  simplified  chemical  kinetics  mechanism,  one  should  pay  particular  attention  to  the 
chemical  reactions  which  maintain  the  energetic  balance  of  the  system.  Detail  of  this  work  is 
presented  in  Appendix  II. 

III.3  Stability  and  Numerical  Analysis  of  Wake  Flows 

The  linear  instability  of  a  family  of  inviscid,  two-dimensional,  variable-density  shear 
layers  and  wakes  is  investigated.  Vorticity  profiles  corresponding  to  a  monotonically  increasing 
velocity  profile  are  first  examined.  A  larger  family  of  initial  vorticity  distributions  which  model 
the  merger  of  two  unequal  vorticity  layers  of  opposite  sign  is  then  considered.  The  latter  is 
obtained  by  superimposing  on  the  former  a  wake  component,  characterized  by  a  spread,  S,  and  a 
velocity  deficit,  W.  The  initial  density  distribution  resembles  a  temperature  spike  and  is 
described  by  a  thickness,  a,  and  a  temperature  ratio,  Tr.  The  stability  propierties  of  the  layers  are 
interpreted  in  terms  of  a  four-dimensional  parameter  space  (W.S.Tr.o)-  The  non-linear  evolution 
of  the  flowfield  is  illustrated  using  the  transport  element  method. 

Flowfield  stability  exhibits  strong  sensitivity  to  the  details  of  the  density  distribution.  In 
the  absence  of  the  wake  component,  the  stability  properties  of  the  heated  layer  are  divided  into 
three  categories  according  to  the  thickness  of  the  density  profile,  <7,  and  the  vorticity  thickness, 
5w  For  <T  >>  5w,  instability  of  the  Kelvin-Helmholtz  mode  in  a  uniform-density  flow  is 
recovered.  When  a~  <v.  die  shear  layer  mode  is  inhibited;  while  this  trend  persists  for  a<  5^, 
the  layer  becomes  characterized  by  the  appearance  of  additional  short-wavelength  unstable 
modes  which  become  dominant  as  a  decreases  and  Tr  increases.  Addition  of  a  wake  component 
is  shown  to  alter  this  behavior,  and  to  oppose  the  stabilizing  effects  of  heat  release.  In  this  case, 
the  shear  layer  mode  always  dominates  the  wake  mode,  and  the  presence  of  heated  sublayer  has  a 


weak  effect  on  the  instability  of  the  vorticity  layer  when  5  is  large,  but  may  influence  the  phase 
speed  of  unstable  waves  whenever  the  zones  of  high  vorticity  and  high  density  gradient  coincide. 
Detail  of  this  work  is  presented  in  Appendix  in. 
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VOPTEX-TRANSPORT  ELEMENT  SIMULATION  OF  THE 
EXOTHERMICALLY  REACTING  SPACIALLY  DEVELOPING  SHEAR 

LAYER 

Marios  C.  Soteriou  and  Ahmed  F.  Ghoniem 

Department  of  Mechanical  Engineering 
Massachusetts  Institute  of  Technology 
Cambridge,  MA  02139 


ABSTRACT 

A  numerical  methodology  is  introduced  which  enables  the  study  of  a  post-transitional 
spacially  developing  exothermically  reacting  shear  layer  over  a  substantial  range  of  the 
governing  parameters.  The  Vortex-Transport  Element  Method,  commonly  used  in  the  simulation 
of  non-reacting  flows  is  further  developed  to  accommodate  an  exothermically  reacting  flowfield. 
The  scheme  is  Lagrangian,  grid-free  and  adaptive  and  solves  the  unaveraged,  time  dependent  and 
coupled  scalar  transport-reaction  and  Navier-Stokes  equations  respectively,  in  their  scalar- 
gradient  and  voracity  forms.  It  exploits  the  Shvab-Zeldovich  formulation  to  provide  solutions 
for  both  moderately  fast  and  infinitely  fast  reactions.  For  finite  reaction  speeds,  Arrhenius 
kinetics  are  used. 

Numerical  results  arc  used  in  a  preliminary  study  of  the  effects  combustion  cxothcrmicity 
on  both  the  flow  and  scalar  fields.  It  is  seen  that  the  externally  forced  flowfield  is  substantially 
modified  in  the  presence  of  combustion  and  a  reduced  growth  is  experienced  mainly  due  to  an 
alteration  of  the  mechanism  by  which  the  vortical  structures,  which  characterize  the  flow, 
interact  This,  together  with  the  decreased  density  within  the  mixing  region  leads  to  decreased 
efficiency  of  combustion. 

The  simpler,  computationally  less  demanding  infinite  reaction  speed  model  is  found  to  be 
an  effective  model  of  the  detailed  finite  reaction  problem  as  long  as  the  reactions  involved  in  the 
latter  arc  fast  compared  to  the  flow,  the  reaction  zone  is  thin,  and  abrupt  transient  phenomena 
such  as  quenching  are  avoided. 


L  INTRODUCTION 

Post-transitional  exothermically  reacting  shear  layers  are  commonly  present  in  many 
combustion  systems.  The  flow,  which  is  a  manifestation  of  the  growth  of  the  instability  of  the 
shear  region  between  two  reacting  fluids  at  different  velocities  represents  an  important 
mechanism  by  which  reactants  mix  and  bum,  in  such  systems.  Experimental  [e  g.  1)  studies 
indicate  that  the  vortical  structures  which  dominate  non-reacting  shear  layers  persist  in  their 
reacting  counterparts  despite  the  substantial  effects  of  the  combusting  field  on  the  flowfleld. 
These  structures  are  found  to  coincide  with  the  region  where  product  exists,  thus  exemplifying 
their  fundamental  importance  to  the  combustion  process.  Combustion,  in  turn,  strongly 
influences  their  evolution  and  interactions  via  the  release  of  chemical  energy  and  the  resulting 
variable  density  field. 

The  effects  of  the  combustion  heat  release  and  of  the  related  variable  density  field  on  the 
flowfield  have  been  the  subject  of  significant  interest  rccendy.  Even  in  the  absence  of  reaction 
the  presence  of  a  variable  density  field  substantially  alters  the  properties  of  the  flow,  modifying 
the  growth  of  the  mixing  region  the  entrainment  from  the  free  streams  and  the  unsteady  evolution 
of  the  eddies  [2].  Experimental  studies  [3.4.1]  have  indicated  that  in  the  presence  of  an 
exothermically  combusting  field  the  shear  layer  growth  is  reduced,  resulting  in  diminished 
efficiency  of  mixing  and  burning.  This  was  initially  a  rather  surprising  finding  since  combustion 
was  anticipated  to  increase  the  size  of  the  mixing  region  via  volumetric  expansion.  Numerical 
studies  mainly  restricted  to  the  simulation  of  temporally  evolving  reacting  shear  layers  (5-7). 
have,  to  some  extend,  been  able  to  reproduce  this  behavior.  In  this  work  the  applicability  of  the 
idealized,  temporally  evolving  flow  model,  in  simulating  reacting  shear  layers  is  fundamentally 
questioned.  The  temporal  model  aims  at  approximating  the  shear  layer  flow  via  a  Galilean 
space-time  transformation:  a  compuutional  domain  is  selected  which  is  moving  with  the  mean 
flow  velocity  and  which  describes  a  small  section  of  the  flowfield.  This  section  is  defined  by  the 
flow  time  scale  in  such  a  way  as  to  include  one  or  two  vortical  structures.  The  small  size  of  the 
domain,  as  compared  to  that  necessary  in  a  spacial  layer  simulation,  is  the  major  source  of  the 
savings  mentioned  earlier.  The  problem  with  the  use  of  temporal  layers  lies  in  the  fact  that  the 
boundary  conditions  in  the  streamwise  direction  (i.e.  at  the  inlet  and  exit)  of  the  domain  are  not 
explicitly  known.  For  this  reason,  artificial  periodic  boundary  conditions  are  imposed.  As  a 
result  the  actual  flow  cannot  be  reproduced  exactly  from  such  a  calculation.  For  example,  the 
asymmetry  imposed  by  the  velocity  ratio,  which  results  in  the  tendency  of  the  uniform  density 
shear  layer  to  intrude  more  into  the  slower  of  the  two  streams,  cannot  be  captured.  Additionally, 
the  streamwise  asymmetry  in  the  shape  of  the  vortical  structures  is  also  lost.  For  the  reacting 
field,  the  limitations  are  more  devastating.  The  periodic  boundary  conditions  remove  any  spacial 
evolution  of  the  combustion-related  properties  of  material  which  crosses  the  domain  in  the 


strcamwise  direction;  only  temporal  evolution  is  permitted.  On  the  other  hand,  the  same  material 
expenences  significant  flow-related  evolution.  This  is  an  uncommon  situation  since  it  requires 
the  flow  speed  to  be  higher  than  the  combustion  speed,  i.e.  high  Karlovitz  number.  For  a  real 
fuel,  characterized  by  a  high  Damkohler  number,  such  a  condition  would  unavoidably  lead  to 
quenching  (termination  of  the  reaction  due  to  excessive  flow  stretch).  In  the  temporal 
simulations  noted  above,  this  problem  does  not  arise  due  to  an  additional,  and  apparently 
independent,  low  Damkohler  number  assumption.  What  is  evident  from  the  above  discussion 
though,  is  that  the  applicability  of  temporal  simulations  is  in  fact  restricted  to  low  Damkohler 
number  combustion  and  unphysical  solutions  would  result  if  higher  values  of  this  parameter  were 
to  be  used. 

Spacially  evolving  shear  layer  numerical  studies  [8,9],  which  attempt  a  much  more 
physical  description  of  the  flow,  have,  on  the  other  hand,  for  the  most  pan  been  restricted  to 
cases  where  the  effects  of  combustion  on  the  flowfield  are  negligible  (low  combustion  heat 
release).  Such  studies  have  been  used  to  discern  the  strucnire  of  the  reaction  zone  arid  the  effects 
of  the  reaction  speed  (Damkohler  number)  on  the  relative  location  of  the  reaction  zone  with 
respect  to  the  large  structures.  It  was  concluded  that  at  small  Damkohler  number,  the  reaction  is 
most  intense  near  the  center  of  the  large  eddy,  while  as  the  Damkohler  number  increases,  the 
reaction  zone  moves  outwards  towards  the  outer  edges  of  the  eddies.  It  was  also  found  that, 
under  conditions  of  unity  stoichiometry,  a  strong  similarity  exists  between  the  products 
concentration  field  and  the  vorticity  field. 

In  contrast,  in  this  work  a  numerical  methodology  capable  of  dealing  with  significant 
combustion  exothermicity  is  pre  ented  and  implemented  in  the  simulation  of  the  forced,  spacially 
developing  shear  layer. 


>  i 


A  two-dimensional,  post-transitional  reacting  shear  layer  is  considered.  Gravitational 
effects  arc  negligible.  Compressibility  effects  are  permitted  under  the  low  Mach  number 
assumption.  According  to  this  assumption,  the  flow  speed  is  much  smaller  than  the  speed  of 
pressure  wave  propagation.  Hence,  spacial  pressure  variations  reach  equilibrium  rapidly  when 
compared  to  the  flow  timescale  and,  thus,  for  thermodynamic  considerations  they  appear 
negligibly  small.  This  can  easily  be  seen  by  considering  that  the  momentum  equation  for  this 
case,  reduces  to  the  statement  that  spacial  gradients  of  the  thermodynamic  pressure  are  equal  to 
zero.  But  while  pressure  variations  are  small  when  compared  to  the  thermodynamic  pressure, 
they  are  not  negligible  when  compared  to  the  other  forces  governing  fluid  motion.  Thus  by 
rescaling  them  with  the  flow  dynamic  pressure  one  can  ascertain  their  effect  on  the  flowfield. 
This  is  expressed  in  the  more  traditional  momentum  equation  used  herein.  (This  distinction 


between  the  thermodynamic  and  dynamic  pressure  in  a  low  Mach  number  combusting  system 
was  initially  proposed  in  rcf.[51  where  a  detailed  derivation  of  the  approximate  equations  of 
motion  is  presented.) 

While,  under  the  low  Mach  number  assumption,  combustion  has  an  insignificant  effect 
on  the  spacial  thermodynamic  pressure  variation,  it  can  substantially  alter  this  pressure  in  time. 
This  is  certainly  true  in  a  constant  volume  domain  where  the  overall  density  is  constant  and 
substantial  pressure  changes  take  place  as  combustion  heat  is  released  and  the  fluid  temperature 
is  raised.  In  an  open  (infinite  volume)  domain,  on  the  other  hand,  combustion  primarily  alters 
the  fluid  temperature  and  density  resulting  to  approximately  constant  pressure.  In  this  work  the 
flow  is  partially  confined  (see  Section  IV).  Nevertheless,  constant  pressure  combustion  is  still 
assumed  since  in  the  cross-stream  direction  the  confining  walls  arc  substantially  far  from  the 
combusting  region  and  in  the  streamwise  direction  the  domain  is  relatively  short.  The  validity  of 
this  assumption  is  further  assessed  from  the  results  by  comparing  the  streamwise  dynamic 
pressure  change  of  each  fluid  stream  to  the  corresponding  inlet  dynamic  pressure.  Thus,  in  this 
formulation  the  thermodynamic  pressure  is  treated  as  a  consunt  both  in  space  and  time  and 
hence  it  is  absent  from  the  non-dimensionalized  equations.  The  pressure  appearing  in  the 
momentum  equation  is  the  flow  dynamic  pressure. 

Combustion  is  assumed  to  take  place  according  to  a  single  step  reaction  which  consumes 
two  reactants,  one  from  each  stream,  to  yield  a  single  product.  The  chemical  kinetics  arc  of  the 
Arrhenius  type.  The  transport  properties  are  constant.  All  species  arc  assumed  to  behave  as 
perfect  gases  with  equal  molecular  weights,  specific  heats  and  mass  diffusion  coefficients.  The 
Lewis  number  is  equal  to  unity.  Diffusion  effects  are  assumed  small  (high  Reynolds  and  Peclet 
number  flow).  This  allows  higher  order  diffusion  mechanisms  to  be  neglected.  Hence,  thermal 
diffusion  (Soret  effect),  second  order  diffusion  terms  in  the  scalar  equations  arising  from 
products  of  density  and  scalar  gradients,  and  heat  production  due  to  fluid  dynamical  viscous 
dissipation,  are  all  neglected. 

Non  dimensionalization  is  carried  out  using  a  length  scale,  Lo,  a  velocity  scale,  Uo  and 
temperature.  To,  and  density,  Po,  scales.  (The  tildes  indicate  dimensional  variables.  In  what 
follows,  absence  of  the  tilde  denotes  a  non-dimensionalized  variable).  The  actual  values  of  these 
scales  are  specified  in  Section  IV.  The  simultaneous  use  of  density  and  temperature  scales 
allows  for  the  different  scaling  of  the  thermodynamic  pressure  as  compared  to  that  of  the 
dynamic  pressure  which  is  scaled  via  the  density  and  the  velocity  scales. 

Under  the  above  assumptions  the  non-dimensionalized  governing  equations  are: 


Continuity 


^  +  pVu  =  0 
dt  ^ 


I 


Momentum 


d!L».ZP +-Lv^u 

dt  p  Re 


Equation  of  State  p  T  =  1 

Chemical  Reaction  ^  — >  (1-h>  )T]p 

Temperature-Species  ^  ^ 


(2) 

(3) 

(4) 

(5) 


with 


and  where 


j 

1  1  2  '  3  ' 

Sj 

t|  Y,|  Y2'|Yp 

Q 

Qo«|  -1 

((+)  Yi  +  Y2  +  Yp=l) 
w  =  Afp2Yi  Y2exp(-^) 


(6) 

(7) 


u  =  (u,v)  is  the  velocity  vector  in  a  right-handed  Cartesian  coordinate  system  x  =  (x.y)  and  t  is 
the  time.  V  ‘s  ‘he  gradient  operator  and  ^  =  ^  +  u- V  the  Lagrangian- material 

derivative,  p  is  the  dynamic  pressure  and  p  and  T  are  the  fluid  mixture  density  and  temperature 
respectively,  tli  denotes  species  i,  i-1,2  being  the  reacting  species  and  i^p  being  the  product 

Yi  =  —  is  the  mass-fraction  for  species  i  (pi  is  the  species’  partial  density)  and  w  is  the  reaction 
P 

rate.  Re  is  the  Reynolds  number,  Re  =  with  v  being  the  kinematic  viscosity.  Similariy,  Pc 

_  V 

is  the  Pcclct  number  Pc  =  a  is  the  thermal  diffusivity,  which  under  the  unity  Lewis 

a 

number  assumption  is  equal  to  the  mass  diffusivity  D  which,  in  turn,  is  the  same  for  all  species. 

Af  is  the  frequency  factor,  Af  =  .  (})  and  0  are  the  molar  and  mass  stoichiometry  ratios 

_  Uo 

respectively  (t  =  ^  H  being  the  molar  mass  of  species  i).  Qo  is  the  enthalpy  of  reaction, 
M2 

Qo  =  where  Cn  is  the  mixture  specific  heat.  T,  is  the  activation  temDcrature  T,  = 

CpTo  P  RTo 

where  E*  is  the  activation  energy  and  R  the  universal  gas  constant 


5 


I 


For  primarily  numerical  reasons  the  above  equations  are  recast  into  non-primitive 
variable  form.  For  the  equations  governing  fluid  motion  this  is  initiated  via  the  use  of  the 
Helmholtz  decomposition  [10], 

Helmholtz  Dec.  u  =  Uw  +  Up  +  u,  (8) 

which  recognizes  that  the  velocity  can  be  split  into  a  voiticity  induced  solenoidal  component,  Uo>, 
(VxU(o  =  Vxu,  V  uo)  =0)  and  two  irrotational  components;  one  induced  by  volumetric 
expansion,  Uc,  and  the  other  by  the  flow  boundary  conditions  Up  ).  The  vortical  component  is 
obtained  from  the  derinition  of  the  voracity  cok  =  Vxu  (k  is  the  unit  vector  normal  to  the  plain  of 
motion)  and  by  the  use  of  the  streamfunction  (y)  i.e. 

Streamfunction  V  y  =  -  ci)  Um  =  Vx(yk)  (9a,b) 

The  irrotational  components  are  obtained  from  the  relevant  continuity  equation  via  the 
introduction  of  the  concept  of  the  velocity  potential  (^): 

Expansion  V^4»e  =  -  -1-^  u*  =  70*  (10a,b) 

p  dt 

"Potential"  V^0p  =  O  Up  =  V0p  (lla,b) 

The  evolution  of  the  voracity  field  required  in  the  evaluation  of  the  vortical  velocity  component 
is  described  by  the  voracity  equation  established  by  taking  the  curl  of  the  momentum  equation  : 

Voracity  ^ijpk  +  (V  u)  (ok  =  ^  V^cok  (12) 

dt  p2  Re 

Density  related  quantities  in  both  equations  (10a)  and  (12)  are  obtained  from  the  temperature 
field  via  the  equation  of  state. 


To  simplify  the  scalar  transport  equations,  Shvab-Zeldovich  (S-Z)  non-reacting  variables 
(X,  Y)  introduced; 

S-Z  Variables  X  =  Yi-0Y2  Y  =  T-  — Yp  (13a,b) 

1+0 

These  variables  are  constructed  from  combinations  of  the  primitive  reacting  scalars  in  such  a  way 
that  the  equations  governing  their  aansport  (obtained  from  algebraic  manipulations  of  equations 
(S))  are  void  of  reaction  terms,  and  are  hence  much  simpler  to  deal  with:  i.e. 


dt  Pe  ^ 


(14) 
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where  P=  X  or  y. 

Depending  on  the  application,  the  particular  choice  of  these  variables  can  also  result  in  further 
simplifications.  A  good  example  of  this  is  illustrated  by  the  7  variable  used  here.  By  combining  * 

the  temperature  with  the  product  mass-fraction  into  one  variable,  one  can  capitalize  on  the  fact 
that  if  the  initial  conditions  imply  that  these  two  properties  are  initially  directly  related  (i.e. 

Y(x,0)s  1),  then  equation  (14)  above  suggests  that  they  will  remain  related  at  all  later  times.  This, 

in  essence,  reduces  the  equations  to  be  solved  by  one,  since  the  solution  for  y  is  trivially  known  * 

(i.e.  7(x.t)=l ). 

The  main  disadvantage  in  the  use  of  Shvab-Zeldovich  variables  is  that  they  impose  the 
substantial  limitations  on  the  choice  of  the  transport  properties  noted  earlier.  This  is  because  in 
order  to  reduce  the  equations,  the  mass-diffusion  coefficients  had  to  be  assumed  equal  for  all 
species  and  the  Lewis  number  had  to  be  chosen  equal  to  unity. 

In  analogy  to  the  treatment  of  the  flow  equations,  the  Shvab-Zeldovich  scalar  transport 
equations  are  recast  into  gradient  form: 

S-Z Gradient  ^  +  g  Vu-t-gx(  oik)  =  J-V^g  (15) 

dt  Pc 

where  g  =  Vp 

Evidently,  integration  of  the  solutions  of  the  gradient  equations  provides  the  scalar  field  only  * 

within  a  constant  This  constant  is  defined  by  the  boundary  conditions  (pp).  Thus,  in  an 
approach  similar  to  the  Helmholtz  decomposition  of  the  velocity  field,  the  total  Shvab-Zeldovich 
scalar  field  solutions  are  obtained  by  adding  these  two  components. 

9 

S-Z  Scalar  Dec.  P^Pg  +  Pp  (16) 

In  the  limiting  case  of  infinite  reaction  rate  combustion,  where  the  reaction  zone  collapses 

onto  a  line  and  reactant  coexistence  is  prohibited,  the  Shvab-Zeldovich  variable  solutions 
together  with  the  equation  of  the  conservation  of  the  species  mass-fractions  (eq.6)  are  able  to  I 

provide  a  complete  description  of  the  reacting  field,  i.e. 


Infinite  reaction:  X,  ^  0 

Yi  =  X,  ¥2  =  0,  Yp=l-X 

(17a) 

A.^0 

Yi=0,  Yz  Yp  =  l-i-^ 

(17b) 

» 

<t>  1> 

The  temperature  field  is  obtained  using  the  definition  of  7.  If  the  reaction  speed  is  finite,  on  the 
other  hand,  the  Shvab-Zeldovich  solutions  and  the  mass-fraction  conservation  equation  do  not 
contain  all  the  information  necessary  to  construct  the  reacting  field.  At  least  one  reacting  scalar  9 
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must  be  explicitly  obtained.  In  the  formulation  presented  here,  this  scalar  is  chosen  to  be  the 
product  mass  fraction,  i.e 

Product  +(1+(|»)^  (18) 

dt  Pc  p 

must  be  solved.  Once  the  product  mass-fraction  solution  is  known,  the  reacting  species  solutions 
can  be  constructed  according  to 


Finite  reaction: 


Yi  = 


Y2  = 


-X  +  d-Yp) 


(19) 


l+<^  l+(^ 

and  the  temperature  solution  by  using  the  definition  of  7,  to  completely  describe  the  reacting 
scalar  field.  It  should  finally  be  noted  that  the  non  reacting  flow  case  is  easily  included  in  the 
finite  reaction  speed  case  by  simply  specifying  Yp  =  0. 


m.  NUMERICAL  SCHEME 

The  numerical  scheme  by  which  the  flow  and  scalar  field  equations  of  the  previous 
section  are  to  be  solved  is  the  Vortex-Transport  Element  Method.  A  non-reacting  version  of  this 
scheme  was  presented  in  ref.  [11]. 

The  numerical  integration  of  the  governing  equations  is  initiated  by  discretizing  the 
vorticity,  "material  density  derivative"  Shvab-Zeldovich  scalar  gradients  and  product 

mass  fraction  via  a  generic  discretization  function.  It  distributes  a  property  ^  over  a  field  of 
elements  which  arc  characterized  by  a  finite  area,  Aj,  and  by  a  strength,  ^i,  locally  distributed  via 
a  radially  symmetric  core  function,  f^.  The  discretization  function  is  also  used  to  reconstruct  the 
discretized  quantities  at  later  times  and  is: 

N 

C(x,t)=2;C.(t)A.(t)f5(lx-Xil)  (20a) 

i=l 

where  Xi  =  Xi(X'0  i*  element  location.  The  core  function  which  is  characterized  by  the  core 
radius  5  within  which  the  most  significant  contribution  of  each  element  is  experienced,  is  a 
second  order  Gaussian,  i.e. 

faW  =  -^xp(-  4^  (20b) 

and  it  enables  a  second  order  accuracy  of  discretization  under  the  condition  that  core  overlap 
between  neighboring  elements  is  maintained  [12]. 
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The  velocity  induced  by  the  discretized  vorticity  field  is  obtained  via  a  discrete, 
desingularized  BiouSavart  law  resulting  from  the  solution  of  equation  (9): 

N 

u<o(*.t)  =  X (21a) 

1=1 

where  r,(t)  =  (i)i(t)A,(t)  is  the  element  circulation  and  K5  is  the  desingularized  kernel  of  the  Biot 
Savart  law  given  by: 

K^.)  =  -|^[l-e«p(-g)l  (21b) 

The  expansion  component  of  velocity  is  obtained  in  a  similar  fashion,  i.e.  via  a  discrete, 
desingularized  equivalent  of  the  convolution  resulting  from  the  governing  Poisson  equation 
(eq.lO): 

N  jq 

i4(x,t)  =  -  S  p  VG5(x-Xi(t))  (22a) 

where  VGj  is  the  desingularized  gradient  of  the  Poisson  equation  and  is  given  by 

VGsd )  =  -  ^  1 1  -  exp(-  ji|2)  (22b) 

The  component  of  the  velocity  field  induced  by  the  boundary  conditions  is  obtained  by 
solving  the  governing  Laplace  equation  (eq.ll)  under  a  Schwartz-Christoffel  conformal 
nuipping  transformation.  The  total  velocity  is  hence  obaincd  via  the  Helmholtz  decomposition. 

The  Shvab-Zeldovich  scalar  solutions  are  obtained  from  the  corresponding  discrete 
gradient  fields  in  an  analogous  manner  since  the  two  quantities  may  also  be  related  via  a  Poisson 
equation,  i.e. 

=  (23) 

Hence  the  Shvab-Zeldovich  scalar  fields  may  be  expressed  as 

N 

P(x,t)  »  X  A.(t)  •  VG5{x-x.(t))  +  Pp  (24) 

i>l 

where  Pp  is  the  integration  constant  obtained  from  boundary  conditions.  The  complete  reacting 
scalar  field  is  obtained  by  using  these  solutions  in  conjunction  with  the  product  mass-fraction 
solutions  which  are  obtained  from  equation  (20). 

The  time  evolution  of  the  flow  and  scalar  fields  is  established  by  numerically  integrating 
the  governing  transport  equations  locally  to  each  element  This  is  done  in  two  fractional  steps. 


'4, 


The  first  step,  the  non-diffusive  step,  includes  all  processes  other  than  diffusion.  During  this  step 
the  element  locations  are  updated  by  numerically  integrating 


that  is,  the  elements  are  advected  with  the  local  velocity  vector.  This  defines  clement 
trajectories,  as  well  as  material  lines.  The  element  vorticity  is  updated  by  locally  integrating  the 
corresponding  circulation  equation. 


dt  “  D 


’‘dt 


were  the  pressure  gradient  has  beea  substituted  by  the  material  acceleration  using  the  relevant 
momentum  equation.  Both  equations  (25)  and  (26)  are  integrated  via  Euler  predictor-corrector 
schemes  and  the  material  acceleration  in  equation  (26)  is  established  by  a  two-step  iteration 
forward-difference  scheme. 

The  direct  integration  of  the  corresponding  Shvab-Zeldovich  scalar  gradient  transport 
equation, 

^+gi  Vu  +  g|x(a),  k)  =  0,  (27) 

is  avoided.  Recognizing  that  under  the  physical  requirements  of  the  problem,  isoscalar  and 
material  lines  may  coincide  and  considering  the  kinemadcal  evolution  of  the  latter  (relating  the 
change  in  length  of  a  material  line  elemental  segment  (Si)  to  the  velocity  difference  at  its  ends. 


^^  =  5l  Vu  ) 

dt 

equation  (27)  may  be  transformed  into  a  simpler  form  [11],  i.e. 


^p8lj 


=  0  with  gi  =  gi  Hi 


(29a,b,c) 


where  H  is  the  unit  vector  normal  to  the  material  line.  Equation  (29)  simply  implies  that  along  a 
material-isoscalar  line 

=  constant  (30) 

p5li 

The  constant  is  specified  by  the  initial  conditions.  51 1  as  well  as  Hi  are  readily  available  due  to 
the  Lagrangian  nature  of  the  scheme  which  trivially  provides  the  topology  and  evolution  of 
material  lines. 

For  the  product  mass-fraction  field  evolution  (finite  reaction  speed  chemistry)  an 
alternative  approach  is  followed.  The  concept  of  product  particles  is  introduced.  Such  particles 
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are  positioned  at  the  center  of  the  computational  elements  and  are  convected  with  the  latter  by 
the  flow.  Product  particles  direcdy  experience  the  combustion  process.  This  is  accomplished  by 
integrating  (via  an  Euler  predictor-corrector  scheme)  the  relevant  transport  equation,  i.e. 


=  (l-Hj))  W 


(31) 


for  each  particle.  The  corresponding  element  strengths  are  then  established  from  the  particle 
values  via  the  discretization  function  (eq.20).  The  effects  of  any  error  associated  with  such  a 
discretization  are  minimized  by  correcting  the  product  particle  values  at  the  end  of  the  complete 
numerical  integration,  i.e.  after  the  diffusion  step,  via  knowledge  of  their  values  at  the  end  of  the 
Hrst  integration  step,  i.e.  after  the  convection-reaction  process. 

In  the  second  fractional  integration  step  diffusion  effects  are  accounted  for.  This  is 
accomplished  via  the  core  expansion  scheme  [13].  Diffusion  is  simulated  by  expanding  the 
element  core  size  according  to: 


5?(t+At)  =  6i^t)  +  4^  (32) 

where  C=Re  for  ^=0)  and  C=Pe  for  ^=g  or  Yp.  This  expression  is  arrived  at  by  using  the  element 
field,  defined  by  the  discretization  function  (eq.(20))  to  analytically  solve  the  governing  diffusion 
equation. 


(33) 

under  *he  assumptions  that  the  strength  and  area  of  the  element  are  constant  and  that  the  core 
radius  is  only  a  function  of  time. 


The  presence  of  the  baroclinic  torque  in  the  evolution  of  the  vorticity  field  requires 
knowledge  of  the  density  gradient  field.  This  in  turn  necessitates  knowledge  of  the  product 
mass-fraction  gradient  field.  In  the  infinite  reaction  speed  case  this  is  trivially  established  from 
the  S-Z  gradients.  For  the  finite  reaction  speed  on  the  other  hand  it  is  established  by 
differentiating  the  product  mass-fraction  field  given  by  equation  (20)  and  is 


where 


N 

VYp(x,t)  =  X  Yp,(t)  Ai(t)  Vf5(x-xi) 

i=l 


(34) 


Vf5(x)  =  -2^fg(x) 
5 


(35) 
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The  severe  stretching  of  the  Lagrangian  mesh  used  in  the  discretization  of  the  transported 
quantities,  which  mcreases  the  distance  between  neighboring  elements,  may  lead  to  the 
deterioration  of  the  accuracy  of  the  solution.  To  overcome  this  problem,  a  scheme  of  local  n«sh 
refinement  is  adopted  whereby  elements  are  continuously  introduced  and  deleted  to  ensure  core 
overlap.  Strong  overlap  is  enforced  near  the  inlet  of  the  domain  by  allowing  a  small  maximum 
separation  between  neighboring  elements.  This  condition  is  relaxed  further  downstream  by 
increasing  this  threshold  value,  thus  allowing  for  efficient  computations  without  compromising 
the  accuracy  of  the  numerical  scheme. 

IV.  FLOW  GEOMETRY.  AND  BOUNDARY  CONDITIONS 

The  geometry  of  the  computational  domain  together  with  some  of  the  boundary 
conditions  are  shown  in  figure  1.  The  shear  layer  evolves  in  a  two  dimensional  channel  of  height 
H  and  length  Xn,»x.  between  two  parallel  streams  (1-top,  2-bottom)  which  mix  downstream  of  a 
thin  splitter  plate.  The  top  and  bottom  walls  are  modeled  as  rigid,  slip,  impermeable  and 
adiabatic  planes  These  conditions  are  satisfied  by  mapping  the  entire  channel  region  on  the 
upper  half  of  a  complex  plane,  and  using  the  appropriate  image  system  of  the  vortex  and 
transport  elements.  At  the  downstream  section,  a  condition  of  vanishing  vorticity  and  scalar 
gradient  is  used  as  outflow  boundary  condition,  and  is  applied  by  removing  the  elements  which 
cross  the  x=Xmax  plane.  At  the  inlet  section,  the  velocity-vorticity,  Shvab-Zeldovich  variables- 
gradients  and  product  mass  fraction  (for  finite  reaction  spwl  simulations)  profiles  are  specified 


as  follows: 

Velocity 

U(x=0,y,t)  =  -t-  crf(^) 

2  2a 

(36) 

Vorticity 

(o(x=0,y,t)  =  -  exp(- 

Vit  o 

(37) 

S-Z  variables 

X(x=0,y,t)  =  ^  +  ^  erf(^^^)  ■)f(x=0,y,t)  =  1 

2  2  O 

(38a,b) 

S-Z  gradient 

V>.(0,y,t)  =  (0 ,  ^0,y,t))  V7(0,y,t)=(0,0) 

(39a,b) 

with 

X  /(y-0.5)2 

— (x=0,y,t)  =  — ^exp(  --) 

By  Vita  o2 

(40) 
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Product 


(41) 


lv-0  5)^ 

.  Yp(x=0.y.t)  = 

In  the  above  profiles,  as  well  as  in  the  numerical  simulations,  the  channel  height  is  used  as  the 
space  non-dimensionalizing  scale,  a  is  the  standard  deviation  of  the  Gaussian  profiles  defining 
the  inlet  voracity,  gradient  of  the  X  variable  and  the  product  mass-fraction.  Its  relation  to  the 
scaling  length  (H)  is  obtained  by  requiring  that  two  wavelengths  of  the  most  unstable  mode  of 
the  uniform  density  shear  layer  (obtained  from  linear  stability  analysis)  fit  within  the  channel 
height.  Finally,  it  is  noted  that  in  the  numerical  simulations  the  velocity  field  is  scaled  with  the 
top  stream  velocity  (UO  whereas  the  density  and  temperature  fields  by  the  common  to  both 
streams  (see  next  paragraph)  values  of  these  properties  (Po.  To ). 

The  profile  for  the  A.  Shvab-Zeldovich  variable  is  obtained  by  assuming  that  each  of  the 
two  fluid  streams  consists  of  a  single  reactant  (i.e.  A.i=l,  X2=-0)  which  at  the  inlet  experiences  an 
errorfunction  type  profile.  Profiles  of  this  type  are  not  unlike  the  experimentally  observed 
profiles  for  the  two  reacting  species.  For  the  y  Shvab-Zeldovich  variable,  the  assumption  is  made 
that  the  temperature  and  product  mass-fraction  profiles  are  directly  related,  yielding  a  constant 
value  proHle.  Under  the  behavior  imposed  by  the  governing  equations  such  an  assumption  is 
valid  for  flows  which  prior  to  combustion  were  of  uniform  temperature  and  in  which  no  product 
was  present. 

In  order  to  avoid  having  to  deal  with  ignition  phenomena  at  the  inlet  for  the  finite 
reaction  speed  case  (note  the  temperature  dependent  nature  of  the  reaction  rate  (cq.7)),  a  finite 
amount  of  product  is  inaoduced  there.  This  is  described  by  the  Gaussian  profile  given  above 
(eq.41)  where  is  chosen  as  0.4.  The  direct  relation  of  the  product  mass-fraction  and  the 
temperature  implies  that  the  temperature  at  the  inlet  is  raised  and  hence  combustion  becomes 
possible.  In  the  infinite  reaction  speed  case  where  reactants  react  on  contact,  ignition  problems 
do  not  exist  and  the  inlet  product  mass-fraction  profile  is  uniquely  defined  by  the  Shvab- 
Zeldovich  variables  as  explained  in  Section  II  (eq.l7).  Hence,  the  scalar  field  inlet  profiles 
specified  above  which  are  used  in  the  numerical  simulation  of  the  flow,  specify  the  inlet  species 
profiles  shown  in  figure  2.  The  related  temperature  and  density  profiles  can  straightforwardly  be 
deduced  fiam  the  product  mass-fraction  profile  as  earlier  explained. 

Initialization  of  the  calculation  is  carried  out  by  assuming  that  the  inlet  conditions  persist 
throughout  the  domain.  Hence,  the  vorticity-S-Z  gradient-product  layer  between  the  two  fluids, 
defined  by  the  above  inlet  profiles  is  discretized  by  distributing  vortex-transport-product 
elements  over  nine  material  layers  (lines)  lying  within  the  support  of  vorticity.  The  elements  are 
of  square  area  of  side  h=0.0195.  The  value  of  the  core  radius  is  6=0.0234  (i.e.5>h). 

Finally,  external  forcing  is  implemented  at  the  inlet.  The  forcing  signal,  shown  in  figure 
3,  consists  of  in-phase  components  of  the  most  unstable  mode  of  the  uniform  density  layer  and 
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its  subhannonic,  both  at  an  amplitude  ApAs=0.025.  The  interaction  of  the  two  forcing 
frequencies  gives  rfsc  to  two  types  of  eddies.  The  eddy  which  forms  during  the  part  of  the  cycle 
which  the  two  components  are  in  phase  -  the  "fundamental"  eddy  -  is  larger  than  the  one  which 
forms  in  the  second  part  of  the  subcycle  -  the  "subharmonic"  eddy  -  during  which  the  two 
components  are  out  of  phase.  The  forcing  is  implemented  by  displacing  elements  at  the  inlet 
according  to  the  forcing  signal. 


♦ 


V.  RESULTS 

Simulations  at  both  finite  and  infinite  reaction  speed  were  carried  out  with  the  double  aim 
of  a)  establishing  a  preliminary  assessment  of  the  effects  of  the  presence  of  combustion  on  both 
the  flow  and  scalar  fields  and  b)  of  providing  a  comparison  between  the  finite  and  infinite 
reaction  speed  models. 

The  time  step  for  all  calculations  was  At=0. 1  and  the  length  of  the  domain  Xmax=5.  The 

IT, 

fluid  dynamical  parameters  were  kept  constant  for  all  runs.  The  inlet  velocity  ratio  is  r=— ^0.5 

and  the  Reynolds  number  (and  Pcclet  number  since  Pc=Rc)  based  on  the  velocity  difference 
across  the  layer  (AU=Ui-U2)  and  on  the  vorticity  layer’s  original  thickness  (d=2o)  is 
Rcd  =  =  5(K).  It  should  be  noted  that  the  Reynolds  number  could  have  also  been  defined 

with  respect  to  the  channel  height  H,  since,  as  was  earlier  explained,  H  is  characteristic  of  the 
largest  cross-stream  scale  of  the  layer.  The  resulting  Reynolds  number  is  Rch  =  =  6400. 

For  the  finite  reaction  speed  cases  the  chemical  paiameters  were  specified  as  follows:  the 

activation  temperature  is  Ta=z^-=  10,  the  frequency  factor  A^=^^^=100  ^=1280), 

RTo  AU  AU 

the  mass  stoichiometry  ratio  0=1  and  the  enthalpy  of  reaction  (referred  to  in  this  text  as  "heat 


diffusion  time  scale  and  tch«=T^xp(^)  the  chemical  reaction  time  scale)  and  the  Karlovitz 

Af  Tf 


number  (Ka=^l**  where  Tnw=-^  is  the  flow  time  scale)  are  Da<j=  1026  and  Kad=  0.49, 
^flw  AU 

respectively  [DaH=  168110,  KaH=  0.038].  The  ranges  of  values  of  these  non-dimensional 
numbers,  which  for  a  given  flow  arc  controlled  by  the  chemical  time  scale,  are  limited  by  the 
available  computational  resources.  This  is  not  only  because  increases  in  the  reaction  speed  will, 
evidently,  necessitate  increases  in  the  temporal  numerical  resolution,  but  also  because  under 
these  conditions  the  reaction  zone  thickness  is  reduced,  thus  requiring  higher  spacial  numerical 
resolution  as  well.  The  values  of  the  Damkohler  and  Karlovitz  numbers  used  here  arc, 
nevertheless,  not  far  from  physically  realistic  values,  describing  a  reasonably  fast  reaction  .  It 
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should  be  pointed  out  though  that  the  same  is  not  necessarily  true  of  the  parameters  used  to 
create  them.  Specifically,  the  values  of  the  activation  temperature  and  the  frequency  factor, 
while  higher  than  those  used  in  previous  work  [7]  are  still  quite  low.  This  is  a  consequence  of 
the  fact  that  due  to  the  nature  of  the  reaction  rate  the  activation  temperature  has  an  independent 
effect  on  the  reaction  zone  thickness.  Even  if  the  effect  of  increasing  the  activation  temperature 
on  the  reaction  time  scale,  is  counterbalanced  by  an  appropriate  frequency  factor,  the  reaction 
zone  thickness  is,  nevertheless,  reduced.  Hence,  spacial  resolution  limits  the  range  of  values  of 
the  activation  temperature  and  the  limitations  on  the  reaction  time  scale,  in  turn,  limit  the 
frequency  factor. 

The  infinite  reaction  speed  model  requires  only  one  chemical  parameter,  the  heat  release. 
This  is  specified  to  be  the  same  as  that  used  in  the  finite  reaction  speed  simulations,  i.e.  Qo=6 

To  assess  the  effect  of  combustion  on  the  flowfield.  the  reacting  calculations  were 
repeated  but  this  time  the  effects  of  the  variable  density  field  on  the  fluid  dynamical  field  were 
ignored.  That  is,  while  the  density  was  allowed  to  vary  due  to  combustion,  the  density  used  in 
the  evolution  of  the  flowfield  was  kept  constant  and  equal  to  its  inlet  value.  Since,  under  the 
earlier  noted  assumptions,  the  variation  in  the  density  field  is  the  only  means  by  which  the 
reacting  field  can  influence  the  flowfield,  then  by  keeping  the  flow  density  invariant  the 
flowfield  remains  ignorant  of  the  presence  of  the  reacting  field.  In  what  follows,  the  cases  where 
the  flow  and  combusting  fields  are  decoupled  in  this  manner,  will  be  termed  as  the  "pnow 
uniform"  cases.  Evidendy,  as  far  as  the  flow  is  concerned  these  cases  ate  indistinguishable  from 
uniform  density  non-reacting  calculations.  Hence,  in  the  analysis  of  the  flowfield  they  will  be 
treated  as  such.  For  the  scalar  field  on  the  other  hand  these  cases  differ  from  uniform  density 
calculations  and  hence  they  will  be  properly  distinguished. 

Figure  4  displays  a  flowfield  comparison  between  the  non-reacting  layer  (left)  and  layers 
reacting  at  finite  (middle)  and  infinite  (right)  reaction  speeds.  The  shear  layer  is  depicted  using 
the  vortex  elements  and  their  local  velocity  vectors.  The  velocity  is  plotted  with  respect  to  the 
inlet  mean  velocity  to  highlight  the  relative  motion  within  the  layer.  To  describe  the  evolution  of 
the  flowfield  in  time,  three  sequential  time  frames  are  shown.  In  order  to  be  able  to  provide  a 
clearer  description  of  the  flow,  which  is  partly  inhibited  by  the  small  size  of  the  flow 
visualizations  of  figure  4,  the  first  time  frame  of  each  case  shown  in  the  figure  is  reproduced  in 
figure  S  at  an  enhanced  scale. 

It  is  seen  that  in  all  cases,  and  in  agreement  with  experimental  evidence,  the  flow  is 
dominated  by  large  scale  coherent  vortical  structures.  The  repetitiveness  of  the  flowfield 
evolution  for  each  case,  a  manifestation  of  the  external  forcing,  is  also  evident  But  the  flow 
behavior  varies  between  the  different  cases  and  comparison  between  them  yields  two  major 


conclusions:  a)  Exothermic  combustion  significantly  modifies  the  flowfield.  resulting  in  reduced 
layer  growth,  and  b)  the  flowfield  characteristics  of  the  finite  and  infinite  reaction  speed  cases, 
while  substantially  different  from  those  of  the  non-reacting  flow,  exhibit  striking  similarities. 

Comparison  betwee:.  the  non-reacting  and  reacting  cases  indicates  that  density  variation 
resulting  from  the  combustion  heat  release  modifies  the  flowfield  by  altering  the  shape,  evoluuon 
and  interactions  of  the  vonical  stnactures.  During  the  initial  rollup,  the  fundamental  eddies 
appear  less  rolled  and  larger  in  overall  size  (area)  than  their  uniform  density,  non-reacting, 
counterparts.  This  increase  in  size  though,  appears  to  be  mainly  in  the  streamwise  direction 
resulting  in  more  elongated,  more  elliptical  structures.  Similar  features  are  experienced  by  the 
subharmonic  eddies  which,  additionally,  appear  much  less  coherent.  Further  downstream,  where 
the  eddies  start  to  interact,  major  differences  are  experienced.  In  the  non-reacting  flow,  each 
subharmonic  eddy  interacts  with  its  downstream  fundamental  eddy.  The  two,  pair  (by  spiraling 
towards  each  other  in  a  clockwise  direction)  to  form  a  larger,  coherent  and  highly  elliptical 
structure  which  continues  to  rotate,  exposing  its  major  axis  to  the  streamwise  flow.  This  results 
in  significant  growth  of  the  shear  layer.  In  the  presence  of  reaction  this  process  is  fundamentally 
altered.  Eddy  pairing  is  resisted  and  the  subharmonic  eddy  appears  to  be  tom  between  its  two 
neighboring  fundamental  eddies  (with  the  downstream  eddy  absorbing  its  larger  part)  in  a  much 
more  continuous  process  than  pairing.  The  resulting  larger  structures  appear  less  elliptic  and 
coherent  than  before  and  they  tend  to  keep  their  major  axis  more  aligned  with  the  streamwise 
direction.  This  impairs  the  layer  cross-stream  growth  significantly.  Finally  it  is  noted  that  the 
convective  speed  of  the  vortical  structures  is  also  altered  in  the  presence  of  reaction.  Figure  5 
displays  this  more  clearly.  Eddies  formed  during  the  same  subcycle  of  the  forcing  function 
(indicated  by  arrows)  appear  to  move  faster  in  the  reacting  flow.  This  feature  is  susp)ected  to  be 
primarily  a  consequence  of  the  volumetric  expansion  associated  with  the  combustion  heat  release 
which  in  semi-confined  flow  like  the  one  considered  here  will  cause  a  streamwise  acceleration. 

The  similarity  of  the  flowfield  between  the  finite  and  infinite  reaction  speed  cases  on  the 
other  hand  is  substantial.  The  eddy  evolution  particularly  at  the  earlier  parts  of  the  domain  is 
strikingly  similar.  Eddy  interactions  are  by  tearing  in  both  cases  and  some  minor  differences  can 
only  be  detected  towards  the  end  of  the  domain  where  for  the  infinite  reaction  speed  case  the 
subharmonic  eddy  appears  to  be  able  to  survive  longer.  The  effect  of  reaction  on  the  eddy  speed 
also  appears  to  be  similar  for  the  two  cases. 

Figures  6  and  7  which  describe  features  of  the  mean  flow  indicate  that  the  above  noted 
instantaneous  behavior  is  repetitive  and  thus  biases  the  mean  flow  characteristics.  Figure  6 
presents  comparisons  of  the  mean  velocity  profiles  (at  a  fixed  downstream  location  (x=3)) 
between  the  non-reacting  flow  and  the  two  reacting  cases;  finite  (6a)  and  infinite  (6b).  The 
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characteristics  of  the  profiles  shown,  are  qualitatively  typical  of  most  downstream  locations 
within  the  computational  domain.  The  steepening  of  the  profiles  in  the  presence  of  reaction, 
which  suggests  smaller  shear  layer  growth,  is  obvious.  The  similarity  between  the  finite  and 
infinite  reaction  case  profiles  is  also  noted.  The  shifting  of  these  profiles  to  higher  speeds,  which 
was  earlier  suggested  by  the  faster  moving  structures,  is  evident.  It  is  interesting  to  note  that  for 
the  reacting  cases  the  mean  velocity  profile  loses  its  monotonic  nature,  typical  of  uniform  density 
shear  layers.  Instead,  close  to  the  fast  free  stream  an  overshoot,  and  close  to  the  slow  free  stream 
an  undershoot  of  the  respective  neighboring  free  stream  values  are  experienced.  This  type  of 
features  of  the  mean  velocity  profile  were  also  documented  in  ref. [6]  for  a  reacting  temporal 
shear  layer  and  were  attributed  to  the  presence  of  the  baroclinic  generation  of  vorticity.  Evidence 
supporting  this  argument  was  also  established  in  our  earlier  work  with  non-reacting  variable 
density  shear  layers  [14],  where  similar  types  of  overshoots  and  undershoots  were  experienced. 

In  figure  7  the  vorticity  thickness  defined  with  respect  to  the  local  free  strearns  as 


5a,= 


Ui(x)  -  U2(x) 


[/3u\ 

L\^yi 

max 

(42) 


(where  the  overbar  indicates  a  time-averaged  property)  and  is,  hence,  representative  of  the 
spacial  cross-stream  thickness  of  the  layer  is  plotted  versus  the  streamwise  coordinate.  (It  should 
be  pointed  out  that  the  definition  of  the  vorticity  thickness  used  here  is  not  necessarily  impaired 
by  the  non-monotonical  nature  of  the  velocity  profile  since  the  overshoots-undershoots  are  small 
(<2%AU).  An  alternative  definition  of  5®  from  the  points  of  the  mean  velocity  profiles  where 
the  velocity  varies  by  5%  of  the  free-stream  value  -i.e.  the  typical  experimental  approach-  would 
give  rise  to  similar  results  since  the  overshoots-undershoots  would  not  be  detected.)  As  in  figure 
6  the  comparison  in  each  part  of  figure  7  is  between  the  non-reacting  and  reacting  flow.  As 
expected,  a  drop  in  the  vorticity  thickness  growth  for  the  reacting  cases  is  clearly  seen.  This  drop 
is  most  significant  towards  the  end  of  the  domain  where  the  effect  of  the  inhibition  of  eddy 
pairing  is  most  pronounced.  Thus,  the  alteration  of  the  eddy  interactions  represents  the  most 
important  mechanism  by  which  the  combustion  exothermicity  reduces  the  forced  shear  layer 
growth. 

Figure  8  displays  the  instantaneous  vorticity  fields  for  the  three  cases  of  figure  5.  It 
should  be  noted  that  only  the  0<x<3  part  of  the  computational  domain  is  indicated.  As  expected, 
the  correlation  with  the  vortical  structures  is  evident.  The  same  is  also  true  for  the  similarity  of 
the  finite  and  infinite  reaction  speed  cases.  But  the  figure  makes  clear  another  difference 
between  the  non-reacting  and  reacting  cases.  The  former  is  characterized  by  vorticity  of  a  single 
sign  (negative)  while  the  latter  by  vorticity  of  both  signs;  a  field  of  mainly  negative  vorticity  with 
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islands  of  positive  vorticity  surrounding  the  vortical  structures.  The  origin  of  the  uniformly 

negative  vorticity  of  the  non-reacting  case  and  the  overwhelmingly  negative  field  of  the  reacting 

cases  should  come  as  no  surprise  since  the  inlet  vorticity  is  negative  (top  stream  is  fastest.  • 

resulting  in  clockwise  rotation).  The  positive  vorticity  of  the  reacting  cases  must  be  a 

manifestation  of  the  baroclinic  generation  since  ft'om  the  mechanisms  of  vorticity  modification 

present  here,  only  this  one  may  change  the  vorticity  sign.  The  presence  of  the  positive  vorticity 

around  the  vortical  structures  must  have  important  effects  on  the  flow  and  is  suspect  for  a  • 

significant  number  of  the  differences  between  the  non-reacting  and  reacting  flowfields.  It 

should,  for  example,  with  its  anticlockwise  rotation,  be  responsible  for  the  overshoots  and 

undershoots  of  the  average  velocity  profiles  as  earlier  suggested.  One  could  speculate  that  its 

location  close  to  the  free  streams  may  inhibit  entrainment  of  iirotadonal  fluid  inside  the  vortical  * 

structures.  By  creating  a  tendency  towards  counterclockwise  rotation  on  the  outskirts  of  the 

eddies  it  may  also  be  responsible  for  the  inhibition  of  pairing,  which  as  was  earlier  noted  is 

initiated  by  a  clockwise  spiraling  of  eddies  around  one  another.  To  be  able  to  clarify  these 

points,  a  closer  scrutiny  of  all  the  mechanisms  by  which  combustion  and  the  resulting  variable  * 

density  field  modify  the  vorticity  field  needs  to  be  carried  out.  This  is  currently  being 

accomplished  and  the  results  will  be  presented  in  a  future  article  (ref.  [15] ). 

The  similarity  of  the  vorticity  fields  of  the  two  reacting  cases,  which  like  earlier  findings 
points  to  the  effectiveness  of  the  infinite  reaction  speed  model  to  reproduce  the  effect  of  the  • 

actual,  finite  reaction  speed  combustion  on  the  fiowfield,  also  provides  an  internal  consistency 
check  of  the  numerical  scheme.  As  noted  in  the  previous  paragraph,  the  dual  sign  of  the  vorticity 
field  of  the  reacting  cases  is  a  manifestation  of  the  baroclinic  generation  of  vorticity.  This 
mechanism  of  vorticity  modification  requires  knowledge  of  the  density  gradient  field.  As  was  * 

earlier  seen  (Section  III)  due  to  the  different  numerical  approaches  in  dealing  with  the  product 
mass-fraction  field  in  the  finite  and  infinite  reaction  speed  cases,  the  density  gradient  is 
established  in  quite  different  ways  for  the  two  cases.  The  qualitative  agreement  of  the  vorticity 
fields  of  figure  8  points  to  the  validity  of  both  approaches.  • 

As  a  final  tool  in  the  investigation  of  the  effect  of  the  combusting  field  on  the  fiowfield, 
statistically  averaged  quantities  are  considered.  In  a  variable  density  field,  such  as  the  one 
resulting  from  the  combustion  heat  release  this  is  usually  done  using  Favre  (density  weighted) 
averaging.  The  advantage  of  this  type  of  averaging  over  the  more  traditional  Reynolds  averaging  i 

is  that  it  results  in  simpler  averaged  equations  by  avoiding  terms  involving  density  fluctuations. 

Its  major  disadvantage  lies  in  the  fact  that  the  resulting  extra  terms  in  the  governing  equations, 
while  fewer  in  number,  arc  less  physically  intuitive.  Furthermore,  they  do  not  incorporate  all  the 
features  of  the  fluctuating  (in  the  Reynolds  sense)  field,  since  fluctuating  quantities  are  part  of  i 


18 


the  Favre  averaged  quantities,  while  at  the  same  time  they  involve  the  effect  of  the  average 
density. 

In  direct  analogy  to  Reynolds  averaging,  the  implementation  of  Favre  averaging  of  the 
equadons  of  modon  results  in  extra  terms  formed  by  the  products  of  the  fluctuadng  quanddes.  In 
the  two  dimensional  case  considered  in  this  work  these  quanddes  are: 

pu  v"  (=pu"v').  pu'\  pv'"*  (43) 

As  pointed  out  in  ref.[16]  the  first  term,  the  turbulent  stress  term,  represents  the  rate  of  transfer  of 
momentum  flux.  The  second  and  third  terms  are  indicadve  of  the  flow’s  turbulent  kinedc  energy 
per  unit  volume  in  the  streamwise  and  cross  stream  direcdons  respectively  and  are  representative 
of  the  turbulence  intensity.  The  total  flow  turbulent  kinedc  energy  can  be  defined  as  : 


KE  = 


pu'^  +  pv"^ 
2 


(44) 


As  noted  earlier,  due  to  the  presence  combustion  heat  release,  the  density  decreases 
significandy  and  this  results  in  decreased  Favre  turbulent  quanddes  (eqs.43  and  44).  It  should  be 
clear  though  that  this  does  not  necessarily  imply  that  the  oscillating  nature  of  the  flow  is 
diminished.  It  only  indicates  to  the  fact  that  the  mean  density  is  reduced.  To  clarify  this  point, 
the  pfiow  uniform  case  with  the  variable  density  of  its  reacting  field  used  in  the  calculation  of  the 
turbulent  quantities,  will  be  provided  for  comparison.  Since  this  latter  case  experiences  the  exact 
same  flowfield  as  the  non-reacting  case  then  the  corresponding  Reynold's  turbulent  quantities  are 
identical.  Hence,  any  difference  in  the  Favre  turbulent  quantities  between  the  two  cases  can  only 
be  a  consequence  of  the  variable  density  and  will  thus  provide  an  indication  of  the  effect  of  the 
decreasing  mean  density. 

Figure  9a  presents  a  comparison  of  turbulent  stress  profiles  between  the  finite  reaction 
speed  case,  its  corresponding  uniform  flow-density  case  (as  explained  above)  and  the  non¬ 
reacting  case.  Figure  9b  presents  an  analogous  comparison  for  the  infinite  reaction  speed  case. 
Profiles  are  displayed  at  a  fixed  downstream  location  (x=2.5)  and  are  typical  of  the  profiles 
experienced  in  the  whole  field.  It  is  clearly  seen  that  the  region  of  significant  fluctuations  is 
related  to  the  shear  layer  region.  This  implies  that  the  shear  layer  represents  the  prime 
mechanism  of  turbulent  mixing  of  the  two  streams.  The  bell-Gaussian  shape  of  the  profiles,  in 
analogy  to  those  of  the  Reynolds  shear  stress  in  a  uniform  density  flow  [17],  clearly  indicates 
that  turbulent  momentum  transfer  from  the  free  streams  to  the  shear  region  represents  the  latter's 
prime  mechanism  of  growth.  But  the  figure  also  clearly  shows  that  this  transfer  is  significantly 
diminished  in  the  presence  of  combustion  heat  release.  Furthermore,  it  is  seen,  that  the  effect  of 
the  decreased  mean  density,  while  significant,  cannot  account  for  the  whole  drop  in  the  turbulent 
stress.  This  suggests  that  heat  release  acts  to  dampen  the  oscillating  nature  of  the  flowfield.  In 
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agreement  with  previous  analysis  the  figure  also  points  to  the  similarity  in  the  characteristics  of 
the  finite  and  infinite  reaction  speed  cases. 

A  similar  picture  is  seen  in  the  turbulent  kinetic  energy  profiles  (figure  10a,b).  They  are 
plotted  at  the  same  downstream  location  as  those  of  the  turbulent  stress.  Again,  it  is  seen  that  the 
presence  of  combustion,  decreases  the  kinetic  energy  of  the  flow  and  that  only  part  of  this 
decrease  can  be  attributed  to  the  decrease  in  the  mean  density;  the  fluctuating  field  is  dampened 
as  well.  It  is  interesting  to  note  the  change  in  the  shape  of  the  reacting  case  profile.  The  single 
peak  profile  experienced  by  the  non-reacting  flow  is  transformed  into  a  three  peak  profile.  The 
fact  that  the  pnow  uniform  case  profile  experiences  a  singe  peak  profile,  suggests  that  the 
transformation  is  not  due  to  the  mean  density  or  its  fluctuation.  Rather  it  is  a  result  of  the  altered 
velocity  field  and  it  implies,  as  is  the  case,  significant  differences  in  the  field's  structure.  The 
shift  from  a  single  to  a  three  peak  profile  invites  the  speculation  that  this  might  be  related  to  the 
alteration  of  the  vorticity  field  in  the  cross-stream  direction  from  an  all  negative  (for  non¬ 
reacting)  to  a  positive-negative-positive  (for  reacting)  field  noted  earlier.  As  was  seen  this 
modifies  the  monotonic  nature  of  the  non-reacting  flow  mean  velocity  profile  to  one  with  a 
gradient  which  experiences  three  peaks  instead  of  one.  Thus,  it  appears  that  the  turbulent  kinetic 
energy  could  be  related  to  the  cross  stream  gradient  of  the  mean  velocity  profile. 

Finally,  a  brief  assessment  of  the  feedback  of  the  modified  flowfield  on  the  reacting  field 
itself  is  to  be  attempted.  Figure  11  displays  two  dimensional  visualizations  of  product  mass- 
fraction  fields.  It  should  be  remembered  that  under  the  assumptions  imposed  on  the  governing 
equations  these  two-dimensional  mass-fraction  maps  are  equivalent  to  temperature  maps  or 
inverse  density  maps.  The  comparison  is  between  the  pnow  uniform  finite  reaction  case  and  the 
corresponding  finite  and  infinite  reaction  spieed  cases.  The  figure  clearly  shows  that  product 
formation  is  strongly  dependent  on  the  evolution  of  the  vortical  structures  which  in  turn  are 
dependent  on  the  destabilization  of  the  vorticity-material  layer  separating  the  two  fluids  at  the 
inlet.  Because  of  this  dependency  figure  1 1  may  also  be  used  to  deduce  the  earlier  noted  effects 
of  the  combusting  field  on  the  flowfield.  The  similarities  between  the  finite  and  infinite  reaction 
speed  cases  are  evident  here  as  well.  This  figure  provides  some  of  the  reasons  why  the  infinite 
reaction  speed  model  is  able  to  provide  such  a  good  approximation  of  the  finite  reaction  speed 
simulations.  It  is  noted  that  the  finite  reaction  speed  result  clearly  indicates  features  of  thin 
reaction  zone  combustion  rather  than  of  distributed  reaction  regions.  Furthermore,  it  is  seen  that 
even  in  the  braids  where  the  strain  is  highest  the  product  concentrations  are  high.  This  implies 
that  quenching,  i.e.  the  collapse  of  the  reaction  zone  due  to  excessive  strained  and  the  associated 
drop  in  temperature,  does  not  take  place.  Both  of  these  effects  point  to  the  fact  that  the 
parameters  chosen  for  the  finite  reaction  speed  simulation  are  such  that  they  result  in  a  fast 


•  • 


20 


reaction  (compared  to  the  flow).  Thus,  it  can  be  concluded  that  the  results  presented  herein 
indicate  that  the  infinite  reaction  speed  model  is  effective  at  mimicking  the  behavior  of  the  finite 
reaction  speed  case  when  the  latter  experiences  a  fast  reaction  (i.e.  approaches  the  infinite 
reaction  speed  limit)  and  not  necessarily  for  alt  other  cases. 

Figure  12  displays  the  quantitative  effect  of  the  combustion  modified  flowfield  on 
product  formation.  This  is  accomplished  by  presenting  the  streamwise  evolution  of  the  product 
thickness  defined  as 

I  <IKx,y)dy 

&t(x)  =  i^ -  (45) 

<t(0,y)dy 

where  <[>  is  the  product  density  pp=pYp.  Hence  the  product  thickness  is  representative  of  the 
mass  of  product  at  a  given  channel  cross-section.  In  each  of  the  two  parts  of  figure  12  the 
product  thickness  of  the  reacting  case  is  contrasted  to  two  versions  of  the  thickness  of  the 
comparison,  pnow  uniform,  case.  It  should  be  remembered  that  the  comparison  case  is  created  by 
simply  decoupling  the  density  of  the  flowfield  from  the  density  of  the  reacting  field  which  is 
allowed  to  vary.  It  should  be  recognized  though,  that  using  a  variable  density  in  the  calculation 
of  the  mass  of  product,  practically  eliminates  the  effect  of  the  change  of  volume  from  this  latter 
quantity.  Thus  the  comparison  between  the  reacting  and  the  pnow  uniform  cases  does  not 
necessarily  provide  a  measure  of  the  relative  difference  in  mass  of  product  but  rather  it  is  more 
indicative  of  the  instantaneous  relative  size  of  the  areas  where  product  exists.  For  this  reason  the 
calculation  of  the  product  thickness  for  the  pnow  uniform  case  is  repeated  using  the  uniform  inlet 
density  for  the  reacting  field  as  well.  This  thickness  is  also  shown  in  figure  12.  Figure  12a 
shows  this  arrangement  of  thicknesses  for  the  finite  reaction  speed  case  and  12b  for  the  infinite 
reaction  speed. 

The  figure  clearly  shows  that  combustion  exothermicity  reduces  product  formation  by 
both  reducing  the  area  where  product  exists  as  well  as  by  reducing  the  density  over  this  area. 
The  second  effect  is  by  far  the  most  dominant.  The  drop  in  the  area  where  product  is  formed  is 
not  as  significant  as  that  of  the  vorticity  thickness  seen  earlier  because  the  latter  only  accounts  for 
the  cross-stream  thickness  of  the  layer,  i.e  the  cross-stream  size  of  the  vortical  structures  but  not 
their  streamwise  size.  As  was  seen  earlier  in  the  flow  visualizations,  for  the  reacting  cases  the 
cross-stream  size  of  the  eddies  is  i  ‘ouced  substantially,  but  their  overall  size  only  marginally. 
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A  numerical  scheme  based  on  the  Vortex-Transport  Element  Method  has  been  presented 
which  is  able  to  subulate  finite  as  well  as  infinite  reaction  speed  combustion  of  significant 
exothermicity  in  a  post-transitional  spacial  shear  layer. 

Numerical  results  indicate  that  the  presence  of  combustion  heat  release  strongly  modifies 
the  flowfield  and,  in  a  forced  layer,  decreases  the  layer  growth  via  an  alteration  of  the  interaction 
of  the  vortical  structures  from  pairing  for  tearing.  The  reduced  growth  in  conjunction  with  the 
reduced  density  within  the  mixing  region  leads  to  diminished  mass  of  products  formed. 

Volumetric  expansion  resulting  from  the  combustion  heat  release  is  seen  to  accelerate  the 
flow  in  the  streamwise  direction.  Baroclinic  vorticity  generation  modifies  the  vorticity  field  in 
such  a  way  that  positive,  counterclockwise  vorticity  appears  at  the  outskins  of  the  vortical 
structures.  The  presence  of  this  positive  vorticity  is  related  to  overshoots  an  undershoots 
experienced  by  the  the  time  averaged  velocity  profiles  as  well  as  to  the  triple  peak  nature  of  the 
Favre-averaged  turbulent  kinetic  energy  profiles.  It  is  suggested  that  a  more  detailed 
understanding  of  the  effects  of  combustion  exothermicity  on  the  flowfield  can  be  obtained  via  a 
closer  scrutiny  of  the  mechanisms  by  which  the  combustion  related  density  field  interacts  with 
the  flow. 

The  infinite  reaction  speed  model  was  found  to  represent  an  effective,  compuudonally 
cost  efficient  approach  in  analyzing  the  effects  of  the  combustion  heat  release  on  the  flowfield  as 
long  as  the  reactions  to  be  replaced  by  the  model  are  fast  compared  to  the  flow. 
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Xmax 


Figure  1.  The  geometry  of  the  computational  domain  together  with  initial  element 
configuration  and  some  of  the  boundary  conditions 


Figure  3.  The  forcing  function  used  to  destabilize  the  flow  and  to  promote  eddy  interactions. 
It  represents  the  cross-stream  distance  by  which  the  elements  are  displaced  at  the  inlet 
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Figure  5.  Flowficld  comparison  (voncx  clement  plots)  for  non-reacting  (top)  and  reacting  at  finite  (middle)  and 
infinite  (bottom)  reaction  speed  shear  layers.  The  time  is  t=10. 


Figure  8.  Vorticity  contours  comparison  between  the  non-reacting  (top)  and  reacting  at 
finite  (middle)  and  infinite  (bottom)  reaction  speed  shear  layers.  The  time  is 
t=10.  0<x<3  pan  of  domain  displayed. 
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Rgure  9.  Comparison  of  the  Favre  turbulent  stress  profiles  between  the  non-reacting, 
reactmg,  and  reacting  with  uniform  density  imposed  on  the  flow  but  not  in  the 
calculation  of  the  stress,  (a)  is  the  finite  and  (b)  the  infinite  reaction  speed  case. 
Streamwise  location  x=2.5 
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Figure  10.  Comparison  of  the  Favre  turbulent  kinetic  energy  profiles  between  the  non¬ 
reacting,  reacting,  and  reacting  with  uniform  density  imposed  on  the  flow  but 
not  in  the  calculation  of  the  stress,  (a)  is  the  finite  and  (b)  the  infinite  reaction 
speed  case.  Streamwise  location  x=2.5 
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equally  subdivide  the  plotting  scale  of  0<Yp<l. 
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Product  thickness  profiles 
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Figure  12.  Comparison  of  the  product  thickness  profiles  between  the  reacting,  reacting 
with  uniform  density  imposed  only  on  fb3w  and  reacting  with  uniform  density 
cases  for  (a)  finite  reaction  and  (b)  infinite  reaction  spe^s.  # 
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Dynamics  of  Unsteady  Strained  Flame 

Abstract 

Transient  response  of  flame  to  an  unsteady  strain  rate  is  analyzed  using  a 
series  of  mathematical  transformations.  Initially,  the  flame  is  located  in  the 
stagnation  plane.  At  time  t=0  ^he  unsteady  strain  rate  is  applied.  Two  basic 
patterns  of  the  strain  rate  are  considered:  a  step  function,  and  a  sinus  wave.  Flame 
response  is  characterized  by  two  parameters:  (1)  burning  velocity  and,  (2)  flame 
location.  The  influence  of  the  Lewis  number  and  the  strain  rate  on  these  two 
parameters  and  on  the  relaxation  time  is  investigated.  For  unity  Lewis  number, 
within  the  range  of  compressive  strains  characteristic  for  a  turbulent  jet  flow  (200 
-  500  1/sec),  the  shapes  of  the  temperature  and  mass  fraction  profiles  remain 
almost  unchanged.  This  leads  to  the  burning  velocity  that  only  slightly  varies 
with  the  strain  rate.  In  the  case  of  non-unity  Lewis  number  the  profiles  are 
significantly  altered  even  by  relatively  weak  strain  rates.  This  happens  due  to  the 
interaction  of  unbalanced  heat  and  concentration  diffusion  and  convection.  The 
changes  in  the  profiles  produce  significant  variation  in  the  burning  velocity.  Some 
analytical  results  are  obtained  for  the  relaxation  time  of  flame  as  a  function  of  the 
strain  rate  and  thermo-chemical  parameters.  The  analytical  derivations  are  based  on 
the  application  of  the  integral  approach  in  the  transformed  domain.  In  order  to 
investigate  the  receptivity  characteristics  of  flame,  the  periodic  strain  rate  is 
applied.  In  the  wide  range  of  frequencies  flame  demonstrates  periodic  response. 
The  average  value  of  the  burning  velocity  is  very  close  to  the  burning  velocity  of 
flame  under  the  average  strain;  phase  shift  between  the  burning  velocity  and  strain 
rate  fluctuations  is  approximately  constant  and  equal  to  -1.3;r .  This  pattern  is 
violated  only  when  the  period  of  the  strain  rate  oscillations  is  much  lower  than  the 
diffusion  time  scale.  The  analysis  of  the  flame  response  suggests  that  the  steady 
state  assumption  can  be  used  with  a  reasonable  accuracy  in  a  flamelet  modeling, 
although  the  phase  shift  should  be  taken  into  account.  Extinction  strain  rate  which 
we  define  as  the  strain  rate  when  the  steady  state  flame  location  crosses  the 
stagnation  plane  for  the  first  time,  is  an  exponential  function  of  the  heat  release. 
This  suggests  that,  in  order  to  get  adequate  values  of  the  extinction  strain  using  a 


simplified  chemical  kinetics  mechanism,  one  should  pay  particular  attention  to  the 
chemical  reactions  which  maintain  the  energetic  balance  of  the  system. 

Nomenclature 

k  -  pre-exponential  factor 

Cp  -  specific  heat  of  mixture  at  constant  pressure  in  J/kg  K 
A  -  thermal  conductivity  of  mixture  in  W/  m  K 
p  -  mixture  density  kg/  m^ 
e  -  strain  rate  of  the  flow  in  1/s 

U  -  velocity  of  the  imposed  stagnation  point  flow  in  m/s 

u  -  velocity  induced  by  the  density  drop  at  the  reaction  front 

Yox,  Yn,  Yf  -  mass  fractions  of  oxidizer,  dilutant  and  fuel 

Cox',  C„,  Cf  -  molar  concentrations  of  oxidizer,  dilutant  and  fuel  in  mole/  m^ 

T  -  thermodynamic  temperature  in  K 

X  -  physical  coordinate  normal  to  the  flame  front  in  m 

Q  -  heat  release  per  mole  of  fuel  in  J/  mole 

W  -  reaction  rate  in  mole/  m^  s 

Ru  -  universal  gas  constant 

EJRu  -  activation  temperature,  in  K 

p  •  pressure  in  Pa 

H  -  enthalpy  of  mixture  in  J/  kg 

^ox,  -  molecular  weights  of  oxidizer,  dilutant  and  fuel  in  kg/  mole 

W^mu  -  molecular  weight  of  mixture  in  kg/mole 

m  -  laminar  flame  eigenvalue 

p  -  non-dimensional  heat  release 

Ze  -  Zeldovich  number 

Un  -  laminar  steady  flame  velocity 

6  -  non-dimensional  temperature 

7]  -  position  coordinate  in  the  transformed  domain,  non-dimensional 

Ab,  Ay,  Ar  -  thicknesses  of  the  temperature,  mass  fraction  and  reaction  rate  profiles, 

respectively,  in  the  transformed  domain,  non-dimensional 

o  -  subscript  denoting  values  in  the  cold  mixture  far  ahead  of  the  flame 

b  -  subscript  denoting  values  in  the  products  region  far  behind  the  flame 

*  -  subscript  denoting  characteristic  values 

Introduction 


A  lot  ot  attention  in  the  last  decade  has  been  paid  to  the  investiaation  of 
tlante  located  in  the  stagnation  point  flow.  It  was  due  to  the  relative  simplicity  of 
the  contiguration,  an  opportunity  to  study  flame/  flow  interaction  and  to  compare 
the  numerical  results  with  the  experiment.  This  configuration  is  interesting  not 
only  from  the  fundamental  point  of  view,  but  also  it  might  be  used  in  the 
modeling  of  a  turbulent  combustion  flow.  The  idea  is  to  split  the  flame  front  in  a 
turbulent  flow  into  a  collection  of  small  planar  flamelets.  Each  of  these  flamelets 
propagates  into  the  unbumed  mixture  with  the  velocity  dependent  on  the  local 
parameters  (effective  flow  strain  rate  etc. ).  Combustion  affects  the  flow  field 
mainly  through  the  change  in  density  near  the  reaction  front  and  through  the  heat 
release.  These  characteristics  depend  on  the  amount  of  fuel  consumed  in  the  flame. 
The  effect  of  flow  on  the  flame  manifests  itself  primarily  through  the  increase  of 
flame  area  which  is  characterized  by  a  parameter  called  stretch  (Law  C.  K.  1988). 
This  characteristic  has  two  components:  due  to  the  strain  rate  and  due  to  the  flame 
curvature.  Study  of  flame  in  a  stagnation  point  flow  helps  to  understand  the 
influence  of  the  strain  rate  on  flame  characteristics,  while  the  effect  of  curvature  is 
accounted  for  on  the  global  level  considering  the  flame  front  as  a  collection  of 
small  planar  flamelets. 

Previous  studies  (  Law  C.K,  1988)  demonstrated  that  in  a  unity  Lewis 
number  mixture  and  in  the  range  of  strain  rates  typical  for  turbulent  flows  ( < 
1000  1/sec.)  the  flame  structure  and  burning  velocity  are  almost  unaffected  by  the 
strain  rate  and  conserve  the  values  corresponding  to  the  unstrained  flame.  This 
leads  to  the  suggestion  that  only  the  evolution  of  the  flame  area  should  be  tracked 
in  a  turbulent  combustion  flow  while  the  burning  rate  should  be  kept  constant 
and  equal  to  the  unstrained  flame  burning  rate  (Meneveau  C.  et.  al.  1991). 

In  this  paper  we  are  studying  a  flame  in  non-unity  Lewis  number  mixture. 
Interaction  of  the  unbalanced  heat  and  concentration  fluxes  with  the  strain  rate 
could  have  a  profound  effect  on  the  burning  velocity.  In  many  cases  this 
characteristic  of  flame  can  not  be  assumed  to  be  equal  to  the  unstrained  value.  We 
utilize  numerical  calculations  to  study  the  response  of  flame  to  the  different 
patterns  of  the  strain  rate.  We  are  interested  in  the  range  of  strain  rates  typical  for  a 
turbulent  flow  and  not  considering  rates  much  higher  than  the  partial  extinction 
strain  (Darabiha  et.  al.  1986)  when  the  steady  state  location  of  flame  is  crossing  the 
stagnation  plane  for  the  first  time. 


Another  important  assumption  that  is  usually  made  in  turbulent  combustion 
flow  simulations  is  the  instantaneous  response  of  flame  to  the  strain  rate 
fluctuations.  In  many  papers  it  has  been  pointed  out  that  this  assumption  is  very 
cmde  ( Rutland  1990i  and  a  transient  model  should  be  applied  in  a  turbulent  flow 
simulation.  In  this  paper  we  are  investigating  the  receptivity  of  flame  to  the  strain 
rate  fluctuations.  The  measure  of  receptivity  is  the  ratio  of  the  averaged  burning 
rate  to  the  burning  rate  under  the  averaged  strain.  If  this  ratio  is  equal  or  close  to 
one  in  a  realistic  range  of  strains  then  the  steady  state  assumption  can  be  safely 
used  although  the  correction  on  the  delay  of  the  burning  velocity  flucmations 
should  be  made.  We  found  that  in  a  practically  important  range  of  strain  rate 
frequencies  the  flame  demonstrates  good  receptivity.  Thus,  for  engineering 
purposes,  the  steady  state  assumption  can  be  used. 

In  Chapter  1  steady  and  unsteady  problems  are  formulated  and  the  solution 
procedures  for  both  of  them  are  described  in  detail.  Solution  of  the  steady 
problem  is  used  to  initialize  the  unsteady  problem.  An  approach  by  Zeldovich 
(1985)  is  applied.  Solution  of  the  unsteady  problem  is  based  on  the  method 
described  by  Rutland  (1990).  This  approach  is  extended  to  include  the  non-unity 
Lewis  number  case.  In  Chapter  2  an  integral  method  is  developed  in  order  to 
obtain  the  functional  dependence  of  the  flame  response  time  on  the  thermo¬ 
chemical  and  flow  parameters.  The  idea  of  the  method  is  to  apply  the  regular 
integral  approach  in  the  transformed  domain  where  the  governing  equations  have  a 
simple  fomi  and  to  use  an  assumption  that  the  ratio  of  the  heat  zone  thickness  to 
the  mass  fraction  zone  thickness  is  a  weak  function  of  time.  In  Chapter  3  a 
numerical  procedure  used  to  solve  the  unsteady  problem  is  described  in  detail.  In 
Chapter  4  we  discuss  the  results  of  calculations  obtained  using  a  set  of  Lewis 
numbers  and  typical  thermo-chemical  parameters  .  Chapter  5  contains 
conclusions. 


1.  Formulation 

In  this  paper  a  plane  flame  is  located  in  the  stagnation  point  flow  produced 
by  two  colliding  jets  of  reactants  and  products  ( Figure  1 ).  X-axis  is  directed 


» 

perpendicular  to  the  flame,  v  -axis  -  along  the  flame.  We  assume  that  the 

externally  imposed  stagnation  point  flow  is  inviscid  and  irrotational.  The  , 

distribution  of  the  x-velociiy  component  in  this  flow  is  given  by  the  following 

equation 

U(x)  =  -£(t)x  (1) 

In  this  expression  £  is  the  diagonal  component  of  the  strain  rate  tensor.  . 


» 


» 

Figure  1  Flame  located  in  the  stagnation  point  flow 
Also,  we  assume  that  the  strain  rate  £  is  a  function  of  time.  . 


Initialization 

Initial  temperature  and  deficient  mass  fraction  profiles  are  obtained  by 
solving  the  problem  of  steady  propagating  premixed  flame.  The  flame  propagation 
is  governed  by  the  following  set  of  equations; 

Continuity: 

-f(.pu)  =  0  (2) 

ax 

Enthalpy: 

puM=-f(XM)  +  QW(y,H)  (3) 

ax  ax  ax 

Deficient  mass  fraction: 
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pu^=  W,  W( Y.H)  1 4 , 

dx  dx  dx 

Equation  of  state: 

p  =  pR^  T/\V^  ( 5 ) 

In  this  problem  one-step  chemical  kinetics  mechanism  is  assumed.  Our  numencal 
simulations  using  full-scale  chemical  kinetics  mechanism  demonstrate  that  the 
bulk  details  of  the  flame/  flow  interactions  can  be  described  using  a  simple 
chemistry  mechanism. 

Reaction  rate: 

W  =  k  Cox  Cf  exp(-EJRT)  (6) 

In  this  paper  we  assume  that  fuel  is  deficient.  The  concentration  of  the  oxidizer  is 
almost  constant  and  equal  to  the  initial  concentration: 

C  -C° 

Using  the  relationship  between  the  mass  fraction  of  fuel  and  its  concentration,  we 
can  rewrite  reaction  rate  expression  as  follows: 

W^l^pY,exp(-EJRT)  (7) 

Wf 

Boundary  conditions  for  these  equations  are: 

jc  =  -oo  H  =  Ho,  jc=:  +  oo//  =  /4,  dHIdx  =  0 

;c=-oo  Y=Yf  ^  jc  =  -i-oo  r  =  0  (8) 

Since  two  second  order  ODE  have  five  boundary  conditions,  one  of  the  parameters 
of  the  problem,  namely  mass  flux  per  unit  area  p  u  can  not  be  specified 
independently  and  is  becoming  an  eigenvalue. 

First,  the  continuity  equation  is  solved.  The  result  of  integration  is: 

p  (jc)  u(x)  =  poUn  =  const 

It  means  that  the  mass  flux  per  unit  area  is  constant  along  the  axis  perpendicular  to 
the  flame  front  in  a  steady  propagating  flame.  Introduction  of  the  following  new 
variables  (Zeldovich  et.  al.  1985)  reduces  the  order  of  the  system: 

^  y  =  -^  f  G  (Hb-Ho)  {WM  (9) 

tit-riQ  L  •  j  ax 

(10) 


H,-Hq 


Here  6  is  reduced  temperature,  y  is  non-dimensional  heat  flux,  m  is  non- 
dimensional  mass  flux,  an  eigenvalue.  Subscript  ♦  denotes  characteristic  values  of 
the  parameters.  We  have  used 

^  =  K,  Cp*  =  Cp,o  =  const  W*  =  yy*'  exp{-EJRTb) 

Wf 
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Enthalpy  equation  wntten  in  terms  of  the  new  variables  takes  the  form: 

-  m  \  +  6)  =  0 

'  dO 


Here  ipiB)  is  defined  as; 

(p(e)  =  {  \V}Jc,)l  (  WXJc,]-  =  4  ^  ~Yj  exp(  Z,  ( 6>- 1 )/( 1 +u(  6-1))) 

and  Ze=  -^(T^-To)  is  Zeldovich  number,  a is  non-dimensional  heat 
_RT,-  T, 

release,  Yf=  Y^  Yf  "  is  non-dimensional  mass  fraction  of  the  fuel.  Another 

equation  can  be  obtained  by  multiplying  original  enthalpy  equation  by  Wf,  original 

mass  fraction  equation  -  by  Q  ,  adding  them  together  and  integrating  the  resulting 

equation.  This  operation  results  in  equation 

puiWfiH  -Hb)  +  QY)=  AWfdtL  +  p  OQdL.  (12) 

Cp  dx  dx 


which  can  be  transformed  into 

n  =  v+  L,=  -A_  (13) 

'  L,  de  pCpD 

Here  Le  is  Lewis  number.  Boundary  conditions  are 

(a)0  =  Oy  =  O  Y=l,@e-ly  =  0  Y  =  0  (14) 

Equations  for  reduced  temperature  (11)  and  mass  fraction  (13)  are  solved 
numerically  using  damped  Newton  algorithm  and  a  shooting  procedure.  Runge- 
Kutta  integration  starts  from  the  hot  side  of  the  flame.  The  initial  values  of  jhe 
derivatives  at  this  point  are  obtained  by  approximating  the  functions  y  and  Y  using 
the  first  terms  in  the  corresponding  Taylor  expansions: 

y  =  /:.(l-0)  T=fc(l-0) 

and  substituting  these  expressions  into  the  governing  equations  to  find  the  values 


of  kx  and  k^.: 


dd 


(0=  1)  =  -/:.= 

(0=1)  = 

de 


1  -  ^ 

/l+-4_ 

1 

/  Le  . 

-ki  =  -  ((^i)  ^  +  mki) 


The  Runge-Kutta  integration  proceeds  up  to  the  point  where  non-dimensional 
reaction  rate  function  (p  (0)  is  equal  to  zero  (reduced  ignition  temperature).  At  this 
point  the  numerical  value  of  the  heat  flux  y  is  compared  with  the  analytical 
solution  of  equation  ( 1 1 )  for  y  where  (pid)  is  put  to  zero  and  a  correction  to  the 
eigenvalue  m  is  made.  The  integration  continues  up  to  the  point  where  the 
analytical  and  numerical  solutions  match.  The  equation  for  the  mass  fraction  is 
satisfied  automatically  (Zeldovich  et.  al.  1985).  Once  the  solutions  for  y  (6)  and 


Y[9)  are  obtained,  the  temperature  profile  in  the  real  domain  is  recalculated  usina 
the  definition  of  the  non-dimensional  heat  Hux  \ . 

K,  =  ',  A„j{k  C.^  exp(-EJRTi,)  Cp  pj  i  (pO/(l-p)+l  kiO/y(  0)  (15) 


Solution  of  the  unsteady  problem 

The  reduced  temperature  and  mass  fraction  profiles  obtained  from  the 
solution  of  the  steady  problem  are  used  as  initial  profiles  for  the  solution  of  the 
unsteady  proulem.  Transient  dynamics  of  one-dimensional  flame  is  governed  by 
the  follov^ing  set  of  unsteady  equations; 

Continuity 

^+^  =  0  (16) 
dt  dx 

Mass  Fraction 


Enthalpy 


Also,  the  same  equation  of  state  (5)  and  reaction  rate  expression  (7)  are  used. 


Boundary  conditions: 

@x=-ooY=Yo,T=To,  @x  =  +  oo  Y=0,T=T,  ,  (19) 

Subscript  0  denotes  the  values  of  the  mass  fraction  and  enthalpy  in  the  cold 
mixture. 

In  this  paper  we  adopted  an  approach  of  C,  Rutland  (1990).  Stagnation  point 
flow  U(x)  is  externally  imposed  on  the  flame.  Due  to  combustion,  a  drop  in 
density  occurs  in  the  flame  zone.  This  density  change  induces  a  velocity 
component  u’(x).  This  velocity  perturbation  vanishes  far  away  from  the  flame 
zone  and  is  non-zero  in  its  vicinity.  The  x-component  of  the  total  velocity  can  be 
approximated  as 

u(x)  =  U(x)  +  u'(x). 

Also,  the  following  boundary  layer  approximations  can  be  made  : 
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d  '  d  i  :  d  ■ 

- - i  »  n: - - ■ 

(Jx  ]  ()\  I  :dz  \ 

The  densit>  lirop  in  the  tlanie  occurs  mainly  along  the  x-axis.  Hence 

;.V.vj  >>  v' .  H’' 

Here  v'  and  vi '  are  perturbations  of  the  velocity  field  in  y  and  z  directions, 
respectively.  Using  these  assumptions,  we  can  rewrite  equations  (16),  (17  )  and  (18) 
as  follows: 


dt  dx  dx 


dV 


dr 


p'^  +  p(U+u)  pD^  -WiY,T)Wf 


dt 


dx  dx 


dY 


dx 


dT  d  I  ,dT 


pc,^+pc,(u*u)  A^|  +  e»'(r,r) 

Ot  ox  uX  \  ox 


(20) 

(21) 

(22) 


These  equations  are  simplified  and  translated  into  a  system  of  reaction- 
diffusion  equations  using  two  transformations.  The  main  ideas  behind  the 
transformations  are  described  in  detail  by  Rutland  (1990)  for  the  premixed  flames 
and  by  A.Ghoniem  and  M.  Soteriou  (1992)  for  the  diffusion  flames. 

First,  the  unsteady  equations  are  non-dimensionalized.  In  the  following 
subscript  's'  denotes  some  characteristic  scale  value,  bar  on  the  top  of  a  variable 
means  that  the  variable  is  non-dimensional.  The  non-dimensional  governing 
equations  are: 

^  +  (JL.  U{x)  -»■  ^  -  0 


dt_ 

dY 


X, 


dx  dx 


U,  mu')  ^^^\pD^\-  -^W{Y,  T)Wj 


1  ? 


dY 


dt 


dx  P  dx 


dT  rTT^~\  t,  t,  I  ^ 
~+  u.  (U+u )  —  —  =  -^  — 

dt  dx  Ps  p  dx 


dx 

aC 

O  dx  j 


pYs 


p.p  CpT. 


W{Y,  T) 


The  first  transformation  (r,  x 


t=  t,  x  = 


-> 
exp  ( 


I  £«')*■ 

Jo 


) 


(23) 


translates  the  non-dimensional  equations  into  a  domain  moving  the  applied  strain 
rate:  _ 

+Uh.B(t) 

dt 


i 


» 


» 


» 


» 


» 


I 


» 


9 


^  >,'A 


^  Bin  ~ 

r)X 


=  4  J-fi  4)  ~ 

P  dx 


pD^\-  H'(K,  n  H  , 

dv  /  pK. 


riT  ~'D,  r,  5 ^  I A 

—7  +  /<,  //  «U)  -tt:  =  —  ■_  —  — 

dr  p,  p  dx  \0  dr 


4(A^U_i^ 


VV  (K,  F) 


a/  •^.  dr  -«/  p,  p  dr\Odr/  p,  p  c,L 

The  part  of  the  convective  term  related  to  the  imposed  velocity  Ufx)  is  eliminated 
in  the  (t,  x  )  domain.  The  next  (Howarth)  transformation 
is  used  to  eliminate  the  part  of  the  convective  term  which  contains  velocity  u  . 

To  apply  the  Howarth  transformation,  we  integrate  continuity  equation  with 
respect  to  x  from  -  oo  to  x  and  take  into  account  that  on  the  products  side,  far 
away  from  the  flame  u’  =  0  . 

The  result  of  the  integration  is 


I- 


dx-^  (pu)  ^0  . 

Xs 


Second  (  Howarth)  transformation  (r,  ^  — >  (/,  rf)  is  defined  as: 


/•* 

J  •  ot 


pdx\ 


t=  t 


Application  of  the  second  transformation  and  equation  (24)  translates  mass  fraction 
and  temperature_equations  as  follows: 

A  =  4  B  =(,)  X  3  ^ .  _Lv^(K  f)  W,  (26) 

3^  x}  Prdr]\  Bt]}  pY. 

^  p®+(£__2_vv(F^  70  (27) 

a;  X.^  P>  BriW  ^  BtjI  Pc,T. 

We  can  make  further  simplification,  assuming  that  Mx)  pix)  =  const 
Now  we  can  rewrite  equations  (26-21)  as: 

^  =  4a,B‘(t)A(^+X_^»'(p;  7),  a=  ^  (28) 

3f  xl  Br]\Bri)  Pc,  T.  c,p. 


L.=  ^  (29 

Br]\Le{r])  BtjI  pY.  D 

A 

The  convective  terms  are  absent  in  the  ( r,  77 )  domain. 

Chemical  time  scale  is  used  as  the  time  scale  of  the  problem  and  heat  diffusion 
length  scale  is  utilized  as  the  length  scale  of  the  problem  : 

t.  =  C/  'lW(Yf,TiB  X.  =  fa^  =  VQb/|e  |  ,  Co  =  PoY//^f 
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Initial  mass  fraction  of  fuel  Y '  scales  mass  fraction;  the  difference  between  the 
flame  and  cold'  temperatures  scales  temperature  T,  =  Ti,  -  T)  .  Then  equations 
(28)  and  (29)  are  rewritten  in  terms  of  the  following  reduced  variables: 

r=  -t-.  f=0=  Lil^ 

Yo  T,  -To 

and  have  the  form 

^  =KaB  ^  Ka  =  \e  |  (30) 

3^  ar/lOr//  P  vT(l,l) 

C  (31) 

dt  dT]\Le{7])  hr]]  P  VT(i,  i) 

Here  we  introduced  Karlovitz  number  to  measure  the  strain  rate  and  used  energy 

conservation  equation  Wf  {H„  -Ho)  =  Q  Yf‘  which  can  be  derived  from  equation 

(12).  Boundary  conditions  are;  _ 

@7/  =  0  r=0,  0=  1  ;  @77  =  +  oo  r=l,0  =  O  (32) 

Expression  for  the  non-dimensional  reaction  rate  is; 

^  ^  Y (33) 

p  mi.  1) 

Parameters  Zg  and  jj.  have  been  introduced  earlier  for  the  steady  propagating 
flame  problem. 


2.  An  integral  solution  in  the  transformed  domain 


Before  solving  numerically  equations  (30),  (31)  we  will  analyze  them  using 
an  integral  approach  and  several  assumptions.  The  purpose  of  this  analysis  is  to 
find  the  thickness  of  the  preheat  zone  as  a  function  of  time.  This  result  gives  us 
the  characteristic  time  scale  of  the  flame  response  as  a  fimction  of  Lewis  number 
Le,  Zeldovich  number  Ze ,  Karlovitz  number /ifa  and  heat  release  parameter  //. 

We  assume  that  the  reduced  temperature  9  and  normalized  mass  fraction  Y 
profiles  in  the  transformed  domain  can  be  approximated  as: 

0(77)=  0,  11>A0  r(77)=  1,  71>Ay 

The  profiles  are  plotted  in  Figure  2.  We  are  considering  the  case  when  Le>  1  and 
the  thickness  of  the  reduced  temperature  profile  in  the  transformed  domain  Ag  is 
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ft 


larger  than  that  of  the  mass  fraction  profile.  The  analysis  for  the  case  o't  Le<  I  is 
similar. 


Figure  2.  Reduced  temperature  and  mass  fraction  profiles  as  functions 
of  transformed  coordinate  Tf,  Le>  I 


We  assume  that  reaction  rate  expression  (33)  can  be  approximated  as 

^  W(Y,^  ^  Y  g  (2,  (0-1))  (35) 

P  mA) 

since  the  heat  release  parameter  fj.  is  of  the  order  of  0.8  and  the  reaction  rate  has 
maximum  close  to  the  point  0=1.  Integrating  equations  (30)  and  (31),  we  obtain: 

^  d rj  (36) 

3(  j„  SrjOiji  I  P  vy(i,i) 

/O 

/•+- 


3(  t 


Integration  of  the  reaction  term  yields 


‘  Jo  Jo  ^  vV(l, 


1) 


d  7]  (37) 


p  ivd.i)  p  vKi.i) 

(38) 

3  Z}  2  Ae  2  Aq 

Integration  of  the  unsteady  and  diffusion  terms  in  equations  (36),  (37)  is 
straightforward: 


» 


» 


> 


» 


» 


•  •  •  • 
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.1,  8  <// 


dr]\drij  2Ae 


^dn  =  -  2-^^ 

dt  8  (it 


dr?  \  dr;/  2  Ay 


Combining  all  terms  together,  we  get  two  expressions  from  equations  (36),  (37); 

+  (40) 

o  at  2  Aff  Ay  Aq 

I ^  o  -7^7^  *  ¥• 

a  at  L  Le  Ay  Ay  Aq 

Here  we  denote 

Ay  Aq  3  Ay  Zg  2  Aq  1  Aq 

Equations  (41)  and  (40)  are  ordinary  differential  equations.  In  order  to  solve  them 
analytically  we  are  making  an  assumption  that  A^  Ay  ~  const .  Actually,  the  ratio 
of  two  thicknesses  is  a  weak  function  of  time.  Now  we  can  rewrite  equation  (40) 
as  follows: 

^  =  8/C,B(f)2+  If  (42) 

dt  3 

The  solution  of  this  equation  with  initial  condition  r  =  0,  Aa  (0)  =  is 


Here 


Aa(/)  =  A^e  *  +  {A^  •  As)  e 


A5=  - ^ - ,  A4^l6fl3 

2  Ka- 16  ff3 


Now,  using  transformation  (25),  temperature  distribution  d(ri)  (34)  and  equation 
of  state  (5),  we  can  obtain  an  expression  for  the  thickness  of  the  reduced 
temperature  profile  in  the  real  domain  Siit)  as  a  function  of  time 


e-^(^^dTj/p(Tj)  j  e-^«'(6!tV(l-/i)+l)dr7  =  e-^a'Aa(l+|-^) 

For  the  mass  fraction  profile  thickness  we  have 


Sr  = 


e-f^<^dTi/p(Ji) )  =  Aa  « z  +  -^  ( z  -  3z 2/4  +  z^/  8)  ,  z  =  Ay/  Aa 


(44a) 


(44b) 


Now  we  can  use  an  approximation 
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(45) 


=  >T7 


and  assume  that  temperature  and  mass  fraction  profiles  in  the  transformed  domain 
are  not  ver>  much  different  from  each  other,  i.e. 

::  =  Ay!  Aq  =  1  +  e,  «  1 .  (46) 

Here  we  implicitly  assumed  that  Lewis  number  is  close  to  unity.  Now  we  can 
substitute  expression  (46)  into  equation  (44b)  and,  retaining  the  first  terms  in  the 
Taylor  expansion  series,  obtain; 

8y  =  e  I  ^  o  ^ 

1  o  I'M  ) 

Substitution  of  this  expression  and  equation  (44a)  into  equation  (45)  yields: 

(47) 

Ay  ^  ^ 

Here  we  also  used  Taylor  expansion  for  Lewis  number,  assuming  that  it  is  close  to 
unity. 

In  a  similar  way  the  ratio  of  the  reaction  zone  thickness  to  the  preheat  zone 
thickness  in  the  transformed  domain.  A,/  Ae,  can  be  evaluated.  Denoting  the 
reaction  zone  thickness  as  5-,  we  obtain: 


^  ~jo  ^  )  =  da e  |2  +  —  ( z  •  3z2/4  +  zV  8)|,  z  =Ar/  Aq  «  1 

Here  z  can  be  used  as  small  parameter.  Taking  this  into  account,  we  obtain: 

5,  =  dae-^<i'z  (l +-^|,  z  =d^da  «  1 .  {<■ 


8r  =  da  e  z  1 1  +  —  |,  z  =  d^  da  «  1 .  (44c) 

In  this  equation  only  linear  in  z  terms  were  retained.  On  the  other  hand, 
reaction  and  preheat  zone  are  related  by  the  following  equation  (Law  C.K.  1988): 

4-^.  (49) 

Using  together  equations  (44a),  (44c)  and  (49),  we  obtain 

,  .  3/i 


da  i  +  JL. 


Taking  into  account  equation  (43),  we  rewrite  equation  (44a)  as  follows: 
Si  =e-^a‘Ai  (l4-^ Ase  +  (A,^  -  As)  e  -^)  (l4  A  oi 


8 


8  l-n 


=(A5  +  (4,^-A5)e-(“«'44)')(l-fi-j^)  . 
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In  this  equation  the  time  scale  of  the  thickness  variation  is  located  in  the 
exponential  term  and  is  equal  to 


r  =  [2K,  - 


ZK,. 


ZK., 


9 


^_L 


1-  exp{  ^Ze 


"  -^d  - 


Using  equations  (47)  and  (50)  to  substitute  for  Ar!  Aq  and  Aq!  Ay  ,  we  obtain  the 
time  scale  of  the  flame  thickness  variation 

T=|2A:a-  (l+\(Lt-l){l+^-—))  a-exp(-^pg(^))(\+^pgip)))\  (51) 

[  9Z2  2  8  1-^  2  2  j 


The  right-hand  side  of  this  equation  consists  of  two  terms:  Karlovitz  number  and  a 
combination  of  Lewis  number,  Zeldovich  number  and  heat  release  parameter.  It  is 
easy  to  demonstrate  that  for  typical  chemical  and  thermal  parameters  the  second 
term  is  usually  smaller  than  the  first.  This  suggests  : 

(1)  Lewis  number  has  minor  influence  on  the  flame  thickness  relaxation  time 

(2)  Relaxation  time  is  inversely  proportional  to  the  applied  strain  rate 
(Karlovitz  number). 

This  conclusions  of  the  integral  approach  will  be  checked  by  the  numerical 
solution. 


* 


■/ 


♦ 


3  .  Numerical  Procedure. 

Equations  (30),  (31)  with  boundary  conditions  (32)  are  solved  numerically 
using  Crank-Nicholson  integration  scheme.  The  source  term  is  linearized  in  order 
of  to  make  the  scheme  implicit. 

In  the  following  discussion,  the  upper  index  denotes  the  time  level,  the 

A 

lower  index  is  the  grid  point  number  in  (r,  rj)  domain.  The  computational  grid  is 
uniform.  Discretization  of  equation  (30)  according  to  the  Crank-Nicholson  scheme 
can  be  written  as  follows: 

i  (RHSr'  +  RHSr),  (52) 

At  2 

Ax  is  the  time  step,  RHS  is  the  notation  for  the  right-hand-side  of  equation  (30) 

RHSr  =KaB  \t  ')^  f  -»■  W  (53) 

B7]\diih  'PM  1^(1,  DA 
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The  equation  is  written  for  n-th  time  evel  and  i-th  grid  point.  We  introduce  the 
following  notation  for  the  non-dimensional  reaction  rate; 

,S4, 

VV'(l.l) 

and  use  an  approximation  of  this  term  on  the  (n-»-l )  tune  level: 

=  a"  +  idi2/dz)" dr  =  a"  +  (dn/d0)^ =  A"  +  (dOJde)?  (dr^-d")  (55) 

Substituting  equations  (53)  -  (55)  into  equation  (52).  we  obtain  a  finite  difference 
approximation  of  equation  (30): 


At 


+  a"  +  (do/dff)"  (er^-e,'')  -t- 


KaBHn-- 


(56) 


A") 


)Af 

dri\drili 

Finite-difference  approximation  of  mass  fraction  equation  (31)  can  be  obtained  in  a 
similar  manner. 

Grid  in  the  computational  domain  (f,  r])  remains  uniform.  Every  time  step 
the  thickness  of  flame  is  compared  with  the  distance  from  the  flame  to  the 
upstream  and  downstream  boundaries  of  the  computational  domain.  If  either  of 
these  distances  is  smaller  than  two  flame  thicknesses,  the  computational  domain  is 
regridded  by  doubling  the  Ar}.  This  procedure  is  similar  to  that  used  by  Rutland 
(1990).  Physical  coordinate  Xi  corresponding  to  the  computational  coordinate  T]i  is 
moving  with  time  according  to  equations  (23)  and  (25).  Each  time  step  a  new 
location  of  x,  is  recalculated  in  such  a  way  that  the  origin  of  the  stagnation  point 
flow,  Xi  =  0  ,  remains  motionless.  For  the  second  derivative  of  6  in  the 
computational  domain  the  following  formula  is  used: 


ft"'  -  2*dr  +  a; 


dr}\d-n}i  ~  (Arif 
We  use  equation  (57)  to  rewrite  equation  (56)  in  the  form  suitable  for  the 
application  of  Thomas  tri-diagonal  matrix  inversion  algorithm  : 

ftci  +Ai  ft  -hAi-i  ft-i  =C,  . 

Expressions  for  the  coefficients  Ai*T‘,  A"',  Ai-i"',  and  G"  are: 

Am  =Am=-  B  ••') 

a:'  =  a:.  B  j- .  I  (diVdB)- 

{Ar]f  At  2 

— T + -Sl-j-a*-  ^(dme):Gr 

2  drj  \  drjh  At  2 
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(57) 


(58) 


(59) 


In  this  formulation  we  use  boundary  conditions  (32),  corresponding  to  the  one- 
flame  configuration.  The  derivative  idfUdG)'  is  evaluated  analytically  using 
equation  ( .'3 ). 


4.  Results  and  discussion 
Set  of  parameters 

In  the  numerical  simulation  we  used  a  set  of  parameters  typical  for  one-step 
methane/  air  combustion: 


Activation  temperature  Ta 
Pre-exponential  factor  k 
Heat  release  parameter  p. 

Initial  temperature  Tq 
Mixture  density  po 
Mixture  specific  heat  Cp 
Mixture  thermal  conductivity  X 
Mixture  Lewis  numbers  Le 


2.0e+04  K 
S.Oe-^OS 
0.83 
300  K 

I.  II  kgtm^ 

J. 12e-^03J/kgK 
0.04  Wtm  K 
0.9, 1,1.1 


These  parameters  were  used  also  in  the  solution  of  the  steady  problem. 


Response  of  flame  to  a  stepwise  variation  of  the  strain  rate 


An  important  information  about  the  flame  transient  characteristics  is 
obtained  considering  the  response  of  flame  to  the  stepwise  variation  in  the  strain 
rate.  Initial  profiles  of  temperature  and  concentration  are  provided  by  the  code 
which  simulates  steady  propagating  premixed  flame.  It  is  described  in  detail  in  the 
previous  chapters.  Strain  rate  e  (r)  is  changing  at  time  t  =  0  from  zero  to  some 
constant  value.  The  response  of  flame  is  characterized  by  two  parameters:  (1) 
burning  velocity : 


Stit) 


and.  ( 2)  integral  flame  location  ; 


Wxdx! 


U'  dx 


Several  remarks  can  be  made  about  the  burning  velocity  (61).  In  one 
dimensional  tube  flow  this  integral  is  equal  to  the  consumption  and  propagation 
velocity  of  the  premixed  flame.  In  the  stagnation  point  configuration  the 
interpretation  of  this  integral  is  less  clear.  If  we  assume  that  after  some  time  flame 
has  stabilized  near  some  location  Xf,  then  we  can  write  mass  conservation 
equation  as  follows 

5u(x)  A{x)  =  Sbjs  {£)  A(xj(e))  (63) 

In  this  equation  A(x)  is  the  area  between  two  chosen  streamlines  at  position  x 
where  the  density  is  equal  to  the  initial  density  po,  'ss  ‘  denotes  steady  state  value 
of  integral  (58),  Sa(x)  is  consumption  velocity  with  respect  to  the  unbumed 
gases.  Equation  (63)  can  be  rewritten  as 

5u(x)  =  Sb.isi£)A(xj(e))/A(x)>Si,.ss(£)  (64) 

for  compressive  strains.  In  the  case  of  2D  plane  stagnation  point  flow  streamlines 
are  given  by  equation 

y(x)  =  C(\i/yx  (65) 

The  area  A(x)  is  reduced  to  2*y(x)  and  equation  (64)  can  be  written  as 

•Sa(jt)  =  Sb.ss(£)x/Xj(£)  (u6) 

Thus,  in  the  stagnation  point  configuration  S„  is  not  only  a  function  of  the  strain 
rate  f  but  also  of  location  x.  In  the  hydrodynamic  calculations  the  thickness  of  the 
preheat  zone  of  flame  is  of  the  order  of  5  mm,  which  is  on  the  verge  of  resolution 
of  a  regular  numerical  scheme.  In  this  case  the  flame  front  virtually  collapses  in 
one  infinitely  thin  surface.  To  calculate  the  evolution  of  this  surface  one  should 
find:  (1)  how  much  of  the  mixture  is  consumed  per  unit  area  and,  (2)  how  the  area 
of  this  surface  is  changing  with  time.  In  the  following  we  are  trying  to  address  the 
first  question. 


Response  of  flame  to  the  step  strain  rate  have  been  obtained  for  parameters 
(60)  and  Lewis  numbers  Le  =  0.9. 1.0, 1.1.  In  Figures  3a  and  3b  strained  flame 
burning  velocity  5*  and  flame  location  Xf  are  plotted  as  functions  of  dimensional 
time  for  different  values  of  stepwise  strain  rate  and  unity  Lewis  number.  The  range 
of  strain  rates  is  from  700  to  1600  1/sec.  When  the  strain  rate  is  increasing,  steady 


Figure 


state  flame  location  is  approaching  stagnation  plane  and.  at  some  value  of  strain, 
crosses  it.  This  value  is  another  characteristic  of  tlame.  Rather  loosely,  we  will 
call  this  \alue  of  the  strain  rate  the  "e.xtinction  strain”.  .Acmally,  the  real 
extinction  doesn't  happen  at  this  strain,  but  the  negative  bummg  velocities  which 
occur  after  the  crossing  of  the  stagnation  plane  are  of  little  interest  for  us  here  since 
they  are  absent  in  a  typical  turbulent  flow.  For  unity  Lewis  number,  when  the 
flame  crosses  the  stagnation  plane,  burning  velocity  is  lower  than  the  unstrained 
value  only  by  8  %.  This  happens  when  the  strain  rate  is  equal  to  1,100  1/sec. 
Strains  higher  than  1,100  1/sec  cause  the  flame  to  move  further  into  the  products 
side  until  it  reaches  the  location  -0.2  mm  (see  Figure  3b) .  An  additional  increase 
in  the  strain  rate  makes  the  profiles  steeper  while  the  flame  itself  starts  to  move 
backward. 

In  Figures  4  and  5  the  same  data  is  plotted  (or  Le  =  l.l  and  Le  =  0.9, 
respectively.  Now  the  flame  response  is  substantially  different.  When  Le  is  lower 
than  unity  (  see  Figure  5),  steady  state  burning  velocity  is  actually  higher  than  the 
corresponding  unstrained  value  although  the  flame  itself  is  still  moving  into  the 
stagnation  plane.  The  extinction  strain  for  Le  <  I  case  is  1,400  1/sec.,  i.e.  21  % 
higher  than  the  value  of  the  case.  The  flame  is  "stronger"  and  more 
resistant  to  external  disturbances.  From  the  profile  point  of  view  this  is  explained 
by  the  particular  structure  of  the  preheat  zone  in  the  Le  <  /  case.  The  flow  is 
modifying  the  preheat  zone  in  such  a  way  that  the  maximum  of  the  reaction  rate  W 
increases  while  the  thickness  of  the  W  profile  decreases  only  slightly.  The 
situation  is  reversed  in  all  aspects  when  Le  =  LI  (see  Figure  4  ).  Now  the 
extinction  strain  is  approximately  900  1/sec.  Burning  velocity  decreases  due  to 
the  drop  in  the  maximum  of  the  reaction  rate  W, 

In  a  realistic  turbulent  jet  the  values  of  the  instantaneous  strain  rates  rarely 
exceed  several  hundreds  1/sec.  Thus,  for  example,  for  very  high  strain  rate  of  700 
1/sec.  steady  state  burning  velocity  in  a  Le  =  /  mixture  will  change  by  -  2.4  %,  in 
Le  =0.9  mixture  it  will  increase  by  +15%,  and,  finally,  in  Le=l.l  mixture  it  will 
drop  by  18%  from  the  respective  unstrained  flame  values. 

We  can  conclude  that  even  for  low  values  of  strain  in  the  case  of  non-unity 
Lewis  number  burning  velocity  is  significantly  modified  due  to  the  convection/ 
preferential  diffusion  interaction.  It  is  not  safe  to  assume  that  the  strained  flame 
burning  velocity  is  approximately  equal  to  the  unstrained  value,  as  it  was  in  the 
case  of  Le  =  I.  In  a  flamelet  model,  burning  velocity  has  to  be  considered  as  a 


Flame  location  [mm] 


0  0.002  0.004  0.006  0.00«  0.01  0.012  0.014  0.016 

Time  [sj 

(b)  Flame  location  as  a  function  of  time,  Lesl.l 


Figure  4  Burning  velocity  and  flame  location  as  functions  of  time,  Le=l.l 


0  0.00S  0.01  0.013  0.02  0.023  0.03  0.033 

Time  [s] 

(b)  Flame  locaUon  as  a  function  of  time,  LesO.9 

Figure  S.  Burning  velocity  and  flame  location  as  functions  of  time,  Le=0.9 
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function  ot  Le  and  f.  In  the  next  pan  of  this  chapter  we  vull  trx  to  answer  the 
question  of  whether  thi>  function  should  depend  on  tune. 

riie  receptis  it\  ot  flame  can  be  characterized  b\  a  rela.xation  time.  ITiis  time 
can  be  introduced  in  seseral  different  vvays  but,  of  course,  the  functional 
dependence  of  it  on  the  thermal  and  chemical  parameters  shouldn't  be  influenced 
by  the  way  of  introduction.  Rutland  (1990)  interpolated  burning  velocity  as  the 
following  function  of  time: 

5*  =  Co  +  Cl  exp  i-t/  T)  i 61) 

and  used  the  time  scale  i  as  a  characteristic  of  flame  response.  Figure  3  (  Le  =  I  ) 
demonstrates  that  this  type  of  interpolation  produces  good  result  only  when  the 
strain  rate  is  higher  than  the  partial  extinction  strain  rate,  i.e.  e  >  1 100  ^  ‘ .  These 
high  values  of  strain  rate  were  mostly  used  in  the  Rutland’s  paper.  For  lower 
values  of  strain  (  see  Figure  3)  flame  demonstrates  fluctuating  response  pattern 
typical  for  almost  critically  damped  system.  In  this  case  we  characterize  the  flame 
response  by  the  settling  time,  i.e.  the  time  where  the  amplitude  of  the  burning 
velocity  fluctuations  around  the  steady  state  value  is  becoming  less  than  0.5  %  of 
the  steady  state  value. 

In  Figure  6  the  settling  time  is  plotted  as  a  function  of  strain  rate  in  log-log 
coordinates  for  several  values  of  the  Lewis  number.  The  purpose  here  is  compare 
the  results  of  numerical  simulations  with  analytical  expression  (51)  for  the 
relaxation  time  as  a  function  of  strain  rate  ( or  Karlovitz  number  Ka  ),  Lewis 
number  Lg,  and  other  parameters  obtained  from  the  integral  analysis  performed  in 
the  transformed  domain  It  is  clear  from  the  Figure  that  unity  Lewis  number 

(  Figure  6b)  is  a  special  case  as  far  as  the  relaxation  time  is  concerned.  Close  to  the 
partial  extinction  the  restructuring  of  the  preheat  zone  is  taken  place  and  the  slope 
of  Tb.5  %  =  To.5  %(£)  curve  increases  almost  two  times:  from  -2.4  before  the 
extinction  strain  to  -0.8  after.  The  latter  value  is  close  to  -0.724  reported  by 
Rutland  et.al.  (1990).  In  his  case  I  he  used  mainly  the  strain  rate  values  higher  than 
the  partial  extinction  value.  Hence,  only  one  slope  has  been  reported. 

Equation  (51)  contains  Zeldovich  number  Zg  in  the  denominator.  The 
typical  value  of  this  number  is  8.  The  second  term  in  the  sum  on  the  right-hand 
side  of  (51)  is  small  compared  to  the  used  Karlovitz  numbers  .  This  manifests 
itself  in  the  weak  dependence  of  the  relaxation  time  t  on  the  Lewis  number.  The 
conclusion  is  supported  by  the  data  of  Figures  6  (a),  (c)  and  (d).  In  these  Figures 
the  relaxation  time  is  plotted  as.a  function  of  strain  rate  forLe  =  1.2,  0.9  ,  0.8. 

The  slopes  of  the  curves  are  almost  identical,  although  the  curves  themselves  are 
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slightly  shitting  to  the  right  when  Lewis  number  is  dropping.  A  remark  should  be 
made  here  that  the  relaxation  time  mentioned  in  this  chapter  is  the  time  related  to 
the  burning  velocity  v  ariation  while  in  Chapter  II  the  time  scale  of  the  tlame 
thickness  variation  has  been  calculated.  These  time  scales  do  not  exactly 
correspond  to  each  other  since  not  only  the  reaction  zone  thickness  is  changing 
under  the  strain,  but  also  the  maximum  of  the  reaction  rate  is  usually  varying. 

Response  of  flame  to  the  periodic  strain  rate 

Other  aspects  of  the  flame  transient  behavior  can  be  investigated  studying 
flame  response  to  the  periodic  strain  rate.  In  the  previous  chapter  we  obtained  that 
if  the  strain  rate  is  changing  in  a  stepwise  manner  it  takes  some  time  for  the  flame 
to  adjust  to  the  new  strain.  The  higher  is  the  strain  the  lower  is  the  relaxation  time. 
For  our  set  of  parameters  the  relaxation  time  was  of  the  order  of  10'^  sec.  (see 
Figure  3).  As  far  as  the  period  of  the  strain  rate  fluctuations  is  concerned,  one 
could  guess  that  in  a  turbulent  flow  the  smallest  period  (the  highest  frequency)  of 
velocity  fluctuation  is  associated  with  the  diffusion  time  scale.  Flow  structures 
with  higher  frequencies  will  dissipate  due  to  the  diffusion.  The  diffusion  time 
scale  is  also  critical  for  the  flame  phenomena  since  it  characterizes  the  process  of 
"feeding"  the  flame.  We  used  a  set  of  the  periods  of  strain  oscillations  proportional 
to  the  diffusion  time  scale.  In  the  code  periodic  strain 

e{t)  —  Cay  "t"  Camp  sin{co{t-  tss ))  (68) 

is  applied  after  the  burning  velocity  reaches  its  steady  state.  In  this  expression  Cav 
is  the  average  value  of  strain.  Camp  is  the  amplitude  value  of  the  oscillations,  r„  is 
the  time  it  takes  to  reach  the  steady  state.  For  the  set  of  parameters  (60)  we  choose 
the  average  strain  Cay  to  be  equal  to  700  j  amplitude  of  strain  oscillations  - 
0.5*£av  and  frequency  corresponding  to  the  periods  of  5,  2, 1  and  0.3  diffusion 
times.  In  Figure  7  burning  velocity  of  flame  is  depicted  as  a  function  of  time  for 
this  set  of  parameters.  In  Figure  7a  the  period  of  strain  oscillations  is  0.3  diffusion 
time.  The  diffusion  time  scale  is  approximately  10'^  sec.  Thus  the  period  of 
oscillations  is  only  slightly  higher  than  the  chemical  time  scale  10'^  sec.  Flame 
initially  demonstrates  periodic  response  which  later  becomes  irregular.  The 
amplitude  of  fluctuations  is  very  low  since  flame  relaxation  time  10'^  sec.  is  much 
higher  than  the  period  of  oscillations  0.3*  10'^  sec,  and  when  the  flame  starts  to 
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respond  to  the  flow  tluctuation  the  tluctuation  already  has  changed  direction.  From 
Figure  "b  i  T  =  1  r..,.  i  to  Figure  7c  (  T  =  2  tj,,  i  the  amplitude  of  the  burning 
velocity  oscillations  is  increasing  and  the  oscillations  themselves  are  becomim’ 
more  and  more  periodic  and  regular.  Fluctuations  of  the  strain  rate  w  ith  the  period 
T  =  5  induce  burning  velocity  oscillations  with  the  e.xtremum  values 
corresponding  to  the  extremum  values  of  the  strain  rate.  For  our  set  of  parameters 
the  average  strain  is  700  1/sec.  and  the  amplitude  is  0.5.  It  means  that  the 
extremum  values  of  the  strain  rate  are  1050  and  350  1/sec.  The  steady  state 
burning  velocity,  corresponding,  for  example,  to  the  strain  rate  1050  1/sec.  is 
approximately  equal  to  0.01 18  m/s  (see  Figure  3).  The  minimum  of  the  burning 
velocity  fluctuations  in  Figure  7d  is  also  equal  to  this  value.  It  means  that  now 
flame  has  enough  time  to  respond  fully  to  the  strain  rate  variations.  For  the  peak 
strain  rate  of  1050  1/sec.  flame  relaxation  time  is  0.002  sec.  (see  Figure  3)  which  is 
now  comparable  with  the  period  of  oscillations  5*/^,/  =  5  *  10  ^  sec.  For  the  lower 
frequencies  flame  demonstrates  fairly  good  receptivity  to  the  strain  rate 
oscillations.  The  phase  shift  between  the  burning  velocity  and  strain  rate 
fluctuations  is  constant  and  close  to  - 1.3  ;r. 

Another  way  to  characterize  the  receptivity  of  a  flamelet  is  to  average  the 
instantaneous  burning  velocity  over  the  significant  amount  of  time  T  until  it  does 
not  change  any  more,  and  to  compare  it  with  the  steady  state  value  corresponding 
to  the  average  strain  e^y,  i.e. 

R((0,eay,£amp)^  Sb  {t,  £it))  dt  !  Sb,  ss^Eav)  (69) 

j'ss 

If  the  burning  velocity  oscillations  have  the  average  value  equal  to  the 
corresponding  steady  state  value,  then  the  ratio  R  is  equal  to  one.  In  Figure  8  /?  is 
plotted  as  a  function  of  frequency  for  Eamp  =  0.5  Eav .  The  value  of  R  is  very  close 
to  1 .02  for  a  wide  range  of  frequencies. 


Figure  8.  Flame  response  to  periodic  strain  rate.  Ratio  of  averaged  burning  velocity 
to  the  burning  velocity  at  the  averaged  strain  rate.  Lesl,  amplitude  =  O.S 

It  demonstrates  that  if  one  is  interested  in  the  average  burning  velocity  over  a 
time  period,  he  can  replace  the  averaging  of  the  instantaneous  burning  velocity  by 
the  averaging  of  the  strain  rate  and  using  the  burning  velocity  at  this  average  strain. 
The  regular,  periodic  pattern  of  the  flame  response  to  the  sinusoidal  strain  rate 
oscillations  demonstrates  that  the  flame  relaxation  time  is  a  function  of  the 
magnitude  of  change  of  the  strain  rate,  but  not  of  the  direction  of  this  change, 
i.e.  flame  behavior  doesn't  show  significant  hysteresis. 


Influence  of  Heat  Release  on  the  Extinction  Strain 

In  this  paper  a  simple  one-step  chemistry  mechanism  is  used.  In  reality 
many  chemical  reactions  are  taken  place  at  the  same  time  and  the  important 
question  is  to  establish  which  chemical  reactions  should  be  taken  into  account,  or 
where  equilibrium  or  steady-state  assumptions  can  be  used.  A  step  in  this  direction 
is  to  determine  how  chemical  (pre-exponential  constant,  activation  energy), 
thermo-chemical  (heat  release),  transport  (Lewis  number)  and  flow  (strain  rate) 
parameters  influence  flame  transient  characteristics.  The  influence  of  Lewis 
number  and  strain  rate  was  described  in  the  previous  chapters.  In  this  chapter  we 
are  considering  the  effect  of  the  heat  release  on  the  extinction  strain. 


* 


Heat  Release  parameter  s  (T^b-T_0)/T_b 


Figure  9  Extinction  strain  and  unstrained  burning  velocity  as  functions  of  heat  release  parameter, 
Le  =  1;  interpolation  of  the  extinction  strain  : 

EexilJ)  =  3.8704e-025  *  ejtp(76.l8*  fi) ;  interpolation  of  the  unstrained  burning  velocity  : 

Ship)  =  2.973e-015  *  exp(37.m*  p) 
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Extinction  strain  is  by  our  definition  the  strain  when  the  steady  state  flame 
location  crosses  the  stagnation  plane  for  the  first  time.  For  a  given  value  of  the 
heat  release  parameter  p  ( p  =  {Tb-  To)/  Tb  ),  the  value  of  the  extinction  strain 
was  determined  numerically  by  trial  and  error  approach,  keeping  all  other 
parameters  (60)  the  same.  The  results  of  these  calculations  are  depicted  in  Figure 
9.  In  this  figure  we  also  plotted  the  unstrained  flame  burning  velocity.  Both 
functions  are  strongly  dependent  on  the  heat  release  and  are  very  well  interpolated 
by  exponential  functions.  While  the  exponential  dependence  of  the  unstrained 
flame  burning  velocity  on  the  heat  release  is  clear  from,  say,  equation  (15),  the 
exponential  growth  of  the  extinction  strain  is  somewhat  a  less  expected  result. 

In  a  regular  flame,  thickness  of  the  reaction  zone  is  about  ten  times  smaller 
then  the  thickness  of  the  preheat  zone.  The  flow  "feels"  the  flame  only  through  the 
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heat  release  \v.hich  takes  place  in  the  reaction  zone  and  through  the  accompanied 
temperature  rise  and  density  drop.  The  obtained  exponential  dependence  of  the 
extinction  strain  on  the  heat  release  demonstrates  that  while  simplifying  chemical 
kinetics  mechanism  special  attention  should  be  paid  to  the  chemical  reactions 
which  determine  the  energetic  balance  of  the  system.  For  example,  for 
hydrocarbons  this  means  that  one  shouldn’t  expect  reliable  dimensional  results 
from  the  flame  model  unless  he  will  take  into  account  energetically  important  CO 
oxidation  reactions. 


5.  Conclusions 

A  series  of  mathematical  transformations  was  used  m  simplify  the  equations 
governing  one-dimensional  unsteady  flame/stagnation  point  flow  interaction.  • 

Influence  of  flow  and  thermodynamical  parameters  on  the  transient  response  of 
flame  was  studied.  It  was  found  that: 

(1)  even  for  a  low  value  of  strain  in  the  case  of  non-unity  Lewis  number 

burning  velocity  is  significantly  modified  by  the  convection/  preferential  * 

diffusion  interaction.  Strained  flame  burning  velocity  is  not  even  approximately 
equal  to  the  unstrained  value,  which  was  the  case  when  Le  ~  1. 

(2)  under  the  periodic  strain,  in  the  range  of  frequencies  and  amplitudes, 

corresponding  to  the  typical  turbulence  flow,  the  ratio  of  the  averaged  * 

instantaneous  burning  velocity  to  the  burning  velocity  under  the  averaged  strain  is 

close  to  one.  It  proves  that  in  this  range  of  frequencies  the  flame  possess  a  high 
degree  of  receptivity  to  the  strain  rate  fluctuations.  The  phase  shift  between  the 
strain  rate  and  burning  velocity  fluctuations  was  approximately  -1.3  ;r.  To 
evaluate  the  instantaneous  value  of  the  burning  velocity,  steady  state  assumption 
can  be  used  if  the  characteristic  time  of  the  strain  rate  change  is  greater  or  equal  to 
the  relaxation  time  of  flame  when  it  is  affected  by  constant  strain  rate  with  the 
amplitude  equal  to  the  strain  rate  variation  (see  Figure  6d).  The  correction  should 
be  made  on  the  phase  shift.  Flame  relaxation  time  depends  mainly  on  the 
magnitude  of  the  strain  rate  change,  but  not  on  the  direction  of  this  change. 

(3)  the  fact  that  the  extinction  strain  is  an  exponential  function  of  the  heat 
release  parameter  dictates  the  simplifications  of  the  chemical  kinetics  mechanism 
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in  such  a  way,  that  the  most  imponant  reactions  from  the  energetic  point  of  view 
should  be  preserv  ed. 

(4)  For  non-unity  Lew  is  number  0.5  'x  settling  time  is  inversely  proportional 
to  the  applied  strain  rate  for  all  range  of  Lewis  numbers.  The  graph  of  unity  Lew  is 
number  settling  time  plotted  versus  strain  rate  has  pre-  and  post-extinctional 
branches  with  the  slopes  of  -2.4  and  -0.8  in  log-log  coordinates.  In  this  sense,  in 
terms  of  the  relaxation  time,  unity  Lewis  number  comprises  a  special  case  . 
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ABSTRACT-  The  linear  instability  of  a  family  of  inviscid.  two>dimensional.  vai  iabi<^deiisity  shear 
layers  and  wakes  is  investigated.  Vonicity  profUes  corresponding  to  a  monolonirally  increasing 
velocity  profile  are  first  exaittiiied-  A  larger  fantily  of  initial  vonicity  distributions  wlu«  h  model 
the  merger  of  two  unequal  vorticity  layers  of  opposite  sign  is  then  considered.  The  latter  is  obtained 
by  superimp«)sitig  on  the  former  a  wake  component,  characterized  by  a  spread,  anr(  a  velocity 
deficit.  IV.  The  initial  density  distribution  resembles  a  temperature  spike  and  is  described  by  a 
thickness.  and  a  temperature  ratio.  Tr.  The  stability  properties  of  the  layers  aie  interpreted  in 
terms  of  a  four*dimensional  parameter  space  GV.  5,  Tr  .<? )  The  non-linear  evolution  of  the  Bow  field 
IS  illustrated  using  the  transport  element  method. 

Flowfield  stability  e.'thibils  strong  sensitivity  to  the  details  of  the  den-^ity  disii  ibutiou.  (n  the 
absence  of  the  wake  component,  the  stability  properties  of  the  healed  layer  are  divided  into  tiuee 
categories  according  to  the  tluckness  of  the  density  profile, tr,  and  the  vinlicily  thickness.  For 
<*  >>  Su,,  instability  of  the  Kelvin-Hclmholu  mode  in  a  uniform-density  fiow  is  recovered.  When 
<T  Sit,,  the  shear  layer  mode  is  inhibited;  while  this  trend  persists  for  ^  <  Su  ,  tlie  layer  becomes 
characterizefl  by  the  appearance  of  additional  short-wavelength  unstable  modes  whicii  become 
dominant  as  o  decreases  and  Tr  increases.  Addition  of  a  wake  component  is  shown  to  alter  this 
behavior,  and  to  oppose  the  stabilizing  cffect.s  of  heat  release.  In  this  case,  the  shear  layer  mode 
!  always  dominates  the  wedte  mode,  and  the  presence  of  heated  sublayer  lias  a  weak  effect  on  the 
•  instability  of  the  vorticity  layer  when  S  is  large,  but  may  influence  the  phase  speed  of  unstable 
waves  whenever  the  zones  of  high  vorticity  and  high  density  gra<Uenl  coincide. 
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I.  lutroduction 

The  evolution  of  high  Reyiiohls  number  cheniically-rearling  free  shear  flows  is  gov¬ 
erned  by  a  large  number  of  fundamental  processes.  1  hese  prores.scs  can  be  described  • 

in  terms  of  the  dynamic  effects  of  conrbustion.  which  leads  to  the  establishment  of 
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an  expansion  flowfield  and  zones  of  sliarp  density  variation,  and  iii  terms  of  intrinsic 
instabilities  of  the  underlying  shear  flow  which  shape  the  evolution  of  the  vortic- 
ity  field  In  most  cases  of  practical  interest,  these  effects  are  closely  coupled,  as  the 
vorticity-induced  convective  field  governs  the  molecular  mixing  proces.ses  and  hence 
modulates  local  reaction  rates,  while  the  evolution  of  the  chenucal  reaction  affects 
the  vorticity  field  through  flow  divergence  field  and  baroclinic  vorticity 

The  description  of  such  flows  is  complicated,  due  to  nonlinear  flow-combustion 
interactions  and  due  to  the  large  number  of  parameters  which  govern  the  relevant 
molecular  and  chemical  |>rocess<'s  This  complexity  often  necessitates  the  construc¬ 
tion  of  simplified  fundamental  models  which  facilitate  the  isolation  of  particular 
interaction  modes,  and  reduce  in  the  number  of  governing  parameters 

Linear  stability  theory  has  proven  to  be  an  important  tool  in  theoretical  studies 
of  reacting  and  heterogeneous  shear  flows  [1-1'2].  In  most  of  these  studies,  heat  re¬ 
lease  mechanisms  associated  with  the  mixing  of  initially  separated  reacting  species 
are  modelled  by  corresponding  temperature  and  density  profiles  which  are  imposed 
on  an  otherwise  homogeneous  sheair  flow  Thus,  while  dynamic  effects  associated 
with  the  presence  of  an  expansion  field  and  the  details  of  the  chemical  reaction  pro¬ 
cess  are  omitted,  spatial  density  (and  temperature)  variation  resulting  from  heat 
release  mechanisms  is  retained,  t'sing  this  approach,  the  essential  stability  proper¬ 
ties  of  the  flowfield  have  been  predicted.  In  particular,  linear  stability  studies  have 
shown  that  the  presence  of  two  or  more  zones  of  different  density  significantly  af¬ 
fects  the  devi'lopinent  of  the  flow.  For  instance,  it  htis  long  been  observed  that  a 
non-unity  density  ratio  alters  the  growth  of  shear  flows  and  influences  the  entrain¬ 
ment  induced  by  the  vortical  structures  embedded  therein  [1-4],  even  when  gravity 
effects  are  weak  [5-6J.  In  chemically-reacting  shear  layers,  the  effects  of  di-nsity  vari¬ 
ation  are  equally  pronounced;  stability  results  indicate  that  flowfield  stabilization 
or  destabilization  tiiay  occur,  depending  on  the  details  of  the  density  and  vorticity 
distributions  [7-9],  Furthermore,  linear  stability  results  show  that  the  presence  of 
zones  of  large  density  variation  may  also  affect  the  nature  of  flow  instalulities  (e.g  , 
by  altering  the  bouiularies  separating  absolute  and  convective  instability  modes), 
and  may  result  in  reshaping  the  global  feature®  of  the  flow  [10-12] 

Unfortunately,  the  application  of  linear  stability  results  to  predict  the  behavior 
or  reacting  flows  has  been  complicated  due  to  the  large  sensitivity  of  the  results 
to  the  initial  vorticity  and  density  profiles.  In  this  work,  this  issue  is  tackled  by 
analyzing  the  stability  of  heterogeneous  flows  for  a  wide  range  of  initial  conditions. 
Initial  density  (and  temperature)  profiles  are  assumed  which  model  the  development 
of  nonpremixed  combustion.  In  addition,  initial  vorticity  profiles  correspotidiiig  to 
symmetric  and  tisymmetric  shear  layers  and  wakes  are  considered.  A  large  number 
of  initial  conditions  is  thus  constructed,  and  flowfield  stability  is  examined  in  a 
four-dimeiisional  parameter  space  which  models  a  large  family  of  reacting  layers 
and  wakes. 

The  stability  problem  is  based  on  linearization  of  the  inviscid  heterogeneous 
flow  equations.  The  formulation  of  this  fluid  flow  problem  is  described  in  Section  2, 
and  computed  results  are  discussed  in  Section  3.  Trends  in  the  behavior  of  the  flow- 
fields  are  established  and  further  examined  by  performing  non-linear  simulation  of 
selected  cases  using  the  transport  element  method  (e  g.  [13- 14|)  .Major  conclusions 
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are  given  in  Section  4. 


2.  Foniiulation 


•2  1  THE  STABILITY  PROBLEM 

In  a  riglit-lianded  coordinate  system  (x.j/),  the  initial  flow  field  is  given  in  tonus 
of  the  steady  flow  velocity,  L  (y).  and  the  imposed  density  profile.  g{y)  The  flow  is 
assumed  two-dimensional  and  inviscid.  and  both  the  perfect  gas  and  the  low  Mach 
number  approximations  of  the  governing  equations  are  employed  [1.V16]  The  mean 
density  profile  is  assumed  to  be  the  result  of  heat  deposition  by  the  reaction,  while 
diffusion  effects  are  neglected.  Lnder  these  assumptions,  the  fluid  flow  is  governed 
by  the  momentum  and  continuity  equations,  respectively  written  as: 


_ 

0-^  =  -Vp 

(1) 

Dt 

(2) 

o 

II 

3 

(3) 

where  u  =  («,  e)  is  the  velocity,  t  is  time,  V  =  (d/di,  d/Oy)  is  the  gradient  operator, 
D/ Dt  IS  the  material  derivative,  and  p  is  pressure.  The  goveriimg  equations  may 
be  recast  in  vorticity  form  by  replacing  Eqs  (1)  and  (3)  by: 


Duj 

'oi 


—  X  Vp 
6~ 


(4) 


where  is  the  vorticity.  This  formulation  will  be  later  used  in  the  simulation  of  the 
noiilinear  evolution  of  the  flow. 

The  stability  properties  of  the  variable-density  shear  flow  are  studied  by  de¬ 
termining  the  temporal  behavior  of  small  amplitude  disturbances  imposed  on  the 
initial  steady  flow.  Tlicsf  properties  are  analyzed  in  terms  of  the  disturbance  cross¬ 
stream  velocity  componeiu  iq  which  is  first  written  in  the  form: 


vo(y,t)  =  i  (i/)exp(iQ(T  -  ct)) 


(0) 


whereQ  is  the  normalized  wavenumber,  taken  to  be  real,  and  c  =  r,  -(-  ic,  is  the 
complex  wave  speed.  The  cros.s-slream  component  i  {y)  obeys  the  modified  Ray  leigh 
equation: 


V 


ft 


q'V 

g{U  -  c) 


(6) 


with  boundary  conditions,  v[y  —  -foe)  ~  exp(-oj/)  and  i  (j/  —  -x)  ~  explop). 
In  Eq.  (6)  and  the  following,  primed  quantities  indicate  differentiation  with  respect 
to  y. 

We  aire  interested  in  determining  whether  waves  of  the  form  expressed  by  Eq. 
(5)  are  unstable,  i.e.  whether  their  growth  rate  ac.  >  0.  To  this  end,  the  eigenvalue 
problem  is  solved  using  a  shooting  technique,  in  which  Eq.  (6)  is  integrated  from  one 
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siJi-  of  the  lay*  r  to  the  otlier  using  a  fourth-order  Ruiige-kutla  prediotor-rorrector 
scheme  [IT]  We  keep  iterating  m  the  complex  eigenvalue  space  c  using  a  secant 
algorithm  until  the  absolute  change  in  c  is  reduced  below  JO"'  This  procedure  ® 
yiehls  dispersion  relations  which  relate  the  growth  rate  ac,  and  the  phase  spet'd 
to  the  disturbance  wavenumber  a 

2  2  INITIAL  (  ONDITIONS 

The  formulation  stability  of  the  stability  problem  is  coiiiph  ltd  by  sptcifyiiig  9 
the  sleatly,  parallel  shear  flow  velocity  and  density  profiles  A  lamily  of  shear  layers 
and  density  profiles  is  constructed  to  model  the  physical  processes  shown  schemat 
ically  in  Fig  1  We  consider  cheiiiically-reactiiig  layers  formed  bt  merging  oxidizer 
and  fuel  slreaiiLs  downstream  of  a  thin  sphiter  plate  or  a  blnIT  body  Imiiiediately 
following  the  tip  of  the  plate,  the  velocity  distribution  may  be  modelled  as  the 
superposition  of  two  Blasius  profiles,  starting  from  a  vanishing  streamwise  velocity  • 
and  increasing  to  the  free  stream  velocities  I  1  and  L  2  (Fig  lb)  Further  down¬ 
stream  the  velocitv  profile  approaches  that  of  a  shear  layer  in  which  the  velocity 
increases  nioiiotonically  from  one  side  of  the  splitter  plate  to  the  oilier  (Hg  la) 

In  the  third  flow  configuration,  the  incoming  stre.ams  are  se  parated  by  a  bluff  body 
(Fig  Ir)  In  this  case,  the  velocity  field  is  regarded  as  the  su|'er|iosn ion  of  two 
shear  lavers,  each  resembling  that  shown  in  Fig  la  In  all  rases,  the  heat  deposited  9 
by  the  initial  development  of  the  chemical  reaction  is  modelled  by  a  spike  in  the 
temperature  profile 

In  order  to  model  the  various  flow  coiifiguartions  antiripati-d  in  Fig  1.  we  start 
with  the  exp'  iiientally  fitted  velocity  profile  for  large  distances  downstream  of  the 
splitter  plate. 

L'iy)  =  r„.  +  '  '  ^  '  tanh(^)  (')  * 

where  (  IS  the  mean  flow  velocity  and  is  the  local  vortiniy  llinkness  The 
local  vorticily  ihickiu'ss  is  chosen  as  characteristic  length  scale,  and  the  Vfhiiiiy  is 
normalize"'  •<i:ch  that  the  reduced  expression. 

[■(y.  -  laiddy)  » 

replaces  t';  (7;.  As  suggested  by  Koch  [18],  a  continuous  family  of  velocity  profiles, 
which  approximates  all  of  the  distributions  shown  in  Fig.  1.  ran  be  obt.ained  by 
modifying  Eq  (9)  by  letting; 

f '(y)  =  ( 1  +  IT)  tanhfy  -  -  (L  tanh(  1/ +  I'i) 

where  IV  is  the  wake-deficit,  and  as  the  displacement  of  the  vorticily  layer  In  Eq  ^ 

(10),  U  and  f'  are  restricted  such  that  IV  <  0.  and  6  >  0.  When  II  =  0.  Eq  (8)  is 
recovered  for  ^  =  0.  and  increasing  the  value  of  i*'  >  0  results  in  a  pure  translation 
of  the  tanh  profile.  On  the  other  hand,  when  IV  <  0.  near  wake  velocity  profiles 
are  approximated  for  large  6.  and  the  double  Blasius  profile  is  iiiutated  when  6  is 

small.  I 

The  free  stream  density  and  temperature  are  chosen  as  a  reference  density  and 
temperature  scales.  Accordingly,  the  normalized  initial  density  profile  is  taken  as 


» 


> 


•  99999999 


b  =  1 

H'  =  u  2 

II 

o 

H  =  0.8 

n  =12 

b  =  1 

H  =  0, 10299 

IT  =  0  24709 

If  =  0  55080 

14  =  0  860922 

b  =  1 

U  =  0  0926 1 

IT  =  0.229S8 

»■  =  0  52204 

11  =  0  821  19 

TABLE  1 


a  Gaussian  profile  with  standard  variation  <t,  and  is  expressed  in  tt  riiis  of  the 
temperature  ratio  Tr  as  follows. 

T  I  • 

P(i/)  =  1 - V— t“xp(— (10) 

/  r  <?’■' 

In  most  shear  flow  applications,  including  shear  lavers,  the  vorticiiv  thickness  is 
usually  larger  than  the  product  thickness,  so  that  it  <  1  However,  we  Jo  not 
enforce  this  restriction  in  order  to  account  for  fluid  flows  c  haraclenzed  bv  higii 
mass  diffusu ities.  or  bluff  body  flows  When  t'  ~  a .  the  flowfielJ  approximates  the 
merger  of  two  layers  of  unequal  density,  a  flow  configuration  that  has  already  been 
analyzed  (see.  eg  [13.19]). 

1  he  temporal  stability  of  the  family  of  variable-density  layers  described  above  is 
investigateil  in  the  four-dimensional  parameter  space  ( IT,  G  Tr .  cr)  We  consider  four 
values  of  the  temperature  ratio,  Tr  —  1,2, -1,  and  8  0,  i  e  we  start  with  a  uniform- 
density  field  and  then  vary  the  temperature  ratio  in  a  range  that  is  representative 
of  most  cln  micaliy -reacting  flows.  For  each  of  these  cases,  the  effects  of  the  wake 
component  (leficir  and  thickness  are  investigated  by  varying  H'  and  b  in  such  a  way 
as  to  approximate  all  the  flow  configurations  shown  in  Fig  1  In  order  to  separate 
the  effects  of  the  strength  of  the  wake  component,  IV,  from  those  associated  with  its 
thickness,  b.  we  alter  the  values  of  H'  and  b  simultaneously  so  that  velocity  profiles 
having  liaving  the  same  maxima  and  minima  are  obtained  for  all  values  of  b.  In 
addition  to  the  lanh  profile  having  U'  =  6  =  0,  we  consider  12  (H',  d)  combinations 
which  are  described  in  Table  I  The  corresponding  velocity  profiles  are  plotted  in 
Fig.  2  for  all  thirteen  cases,  the  shear  layer  profile  being  included  with  the  set  of 
profiles  having  <5  =  1. 

Thus,  li'  varies  in  a  wide  range  of  wake  deficits,  while  increasing  b  from  I  to  3 
repre.sents  a  migration  from  an  a.symmetric  shear  layer  velocity  disi ribntion  to  an 
asymmetric  wake  profile,  in  which  vorticity  layers  are  well  separated  1  he  stability 
properties  of  variable-density  layers  are  first  determined  assuming  eipial  density 
and  vorticity  thicknesses,  cr  =  1  The  effect  of  the  thickness  of  the  di  nsity  (irofile  is 
then  investigated  by  repeating  the  analysis  at  Tr  =  4  0  for  four  additionalir  values, 
s  =  1  5.0  To.  I)  ).  and  0  25 

2  3  BRIF.F  THFORETK  AL  REVIEW 

The  linear  stability  problem  of  inviscid  incompressible  parallel  shi  er  ITav  has 
been  studied  extensively  (e  g  [10,20-21])  In  this  section,  we  summarize  aspects  of 
the  theory  which  directly  affect  our  search  for  unstable  eigenfunctions  We  first  n'  te 
that  the  classical  stability  results  expressed  by  Ray  leigh  s  theorem  [20-21], 
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A  necessary  condition  for  instability  is  that  the  profile  L’(y)  admits  an  inflection 
point 

and  by  Fjortoft  s  extension. 

A  nece.jsary  condition  for  instability  is  that  t"(f'  —  I’,)  <  0  somewhere  the 
flow,  where  i',  is  the  velocity  at  the  inflection  point 

have  to  be  modified  when  considering  a  variable  density  flow.  It  may  easiliy  be 
shown  that  the  appropriate  generalizations  of  the  above  results  may  be  respec¬ 
tively  expressed  as; 

A  necessary  condition  for  instability  is  that  U"  +  (p  lp)U'  admits  a  sign  change 

and 

A  necessary  condition  for  instability  is  that  I’"  -f-  (p'lp)[''[U  -I',)  <  0  some¬ 
where  in  the  flow. 


Thus,  in  a  variable-density  field,  the  behavior  of  the  quantity  L'"  +  {p' /p)i”  replaces 
that  of  U"  in  the  determination  of  the  stability  of  flow.  .As  indicated  in  Refs.  1  and 
8,  this  leads  us  to  expect  strong  interactions  between  the  density  variation  and  the 
shear  flow,  whenever  zones  of  high  vorticity  and  density  grtidieiU  coexist. 

The  behavior  of  the  curves  of  U"  -¥(p' / p)V'  for  the  flow  configurations  of  Table 
I  is  plotted  in  Fig.  3  for  6  =  1,0-  =  1,  and  T,  =  1,2,4,  and  8,  and  in  Fig  4 
for  S  =  l.Tr  =  4.0  and  a  =  1.5,0.75,0.5,  and  0  25.  By  the  preceding,  we  are 
led  to  expect  unstable  modes  whenever  the  curves  intersect  the  zero  axis.  In  the 
uniform-density  case,  a  single  intersection  point  is  observed  for  the  tank  profile, 
and  two  for  the  asynunetric  layers.  This  is  not  surprising,  since  we  expect  to  be 
observe  one  unstable  shear  layer  mode,  and  one  unstable  wake  mode.  However, 
for  large  Tr^  several  intersection  points  appear  for  the  tank  profile,  and  may  yield 
additional  instability  modes.  This  unexpected  result  is  further  investigated  in  the 
following  section  where  the  behavior  of  these  modes,  whose  appearance  depends  on 
the  vorticity-density  configuration,  is  computed. 

VVe  conclude  this  section  by  extending  to  variable-density  flow  the  analysis  of 
Drazin  aiul  Howard  [22]  who  studied  the  the  long  wave  behavior  of  unstable  modes, 
i.e.  the  limiting  behavior  of  unstable  eigenvalues  as  o  —  U  Our  analysis  is  reduced 
to  a  form  similar  to  the  incompressible  case  by  rew-itmg  the  modified  Rayleigh 
equation  as: 

(11) 

whore 

Z^  =  e(U-c?  (12) 

and  (t>  is  the  perturbation  potential.  Noting  that  Eq.  (12)  is  identical  to  Eq  (16) 
of  Drazin  and  Howard  [22],  we  are  able  to  carry  out  a  similar  analysis  to  the  one 
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performed  there  (Details  »■  dl  be  presented  elsewhere).  In  the  limit  a  —  0,  Eq  (14) 
implies  that  Z^F'  is  constant  and  that  this  constant  is  zero  in  order  to  satisfy  the 
boundary  conditions.  Hence  F  is  constant  in  intervals  where  Z  does  not  vanish,  but 
may  have  jumps  when  Z  =  0  Assuming  F  has  no  jumps,  a  siiiular  argun«'nl  to 
that  presented  in  [22]  shows  that  the  limiting  eigenvalues  satisfy: 

^L+2’!^  =  0  (13) 

For  equal  free  stream  densities,  Eq.  ( 14)  yields  the  unstable  eigenvalue  c  —  i  which 
IS  recognized  as  the  limiting  eigenvalue  of  the  shear  layer  mode  Density  variation 
does  not  seem  to  affect  the  tisymptotic  behavior  of  these  long  waves. 

A  similar  result  is  reached  if  F  admits  jumps.  While  the  details  of  the  algebraic- 
manipulations  are  more  involved  than  in  the  uniform-density  case,  we  are  still  able 
to  show  that  if  F  has  a  jump  at  then  f  ’(j/o)  =  0.  Thus,  velocity  maxima  and 
minima  are  expected  to  represent  limiting  values  of  unstable  eigenfunctions  For 
H  <  0,  the  profiles  considered  in  this  study  admit  a  velocity  minimum,  which  is 
recognized  as  the  limiting  eigenvalue  of  the  (unstable)  wake  mode. 

The  long  wave  approximation  estimates  are  used  in  the  following  section  to 
initialize  the  search  for  unstable  eigenfunctions  and  to  characterize  the  additional 
inflection  points  which  appear  for  the  ianh  profile  at  high  temperature  ratio  (Fig. 
3).  Should  these  additional  inflection  points  correspond  to  unstable  modes,  then  the 
above  argument  shows  that  the  associated  instability  mode  affects  shorter  wave¬ 
length  instability,  since  the  asymptotic  behavior  of  long  waves  Ls  solely  dependent 
on  the  details  of  the  velocity  profile. 

3.  Results 

3  1  STAUIEITY  OF  VARIABLE-DENSITY  LAYERS  AND  WAKES 

Stability  analysi.s  of  the  variable-density  parallel  asymmetric  shear  flow  is  first  con¬ 
ducted  for  the  velocity  profiles  of  Table  I,  and  density  profiles  having  =  ),  and 
Tr  =  I  '2. 4,  and  8.  Since  a  sharp  estimate  of  an  upper  bound  on  unstable  wavenum¬ 
bers  was  not  sought,  a  complete  search  for  the  unstable  eigenfunctions  over  a 
bounded  wavenumber-eigenvalue  region  cannot  be  easily  conducted.  Instead,  the 
search  for  unstable  waves  is  initiated  by  concentrating  on  the  behavior  of  long 
waves,  and  extrapolating  the  asymptotic  behavior  of  unstable  shear  layer  and  wake 
modes,  as  determined  by  the  theoretical  predictions  To  this  end,  the  wavenum¬ 
ber  is  increased  incrementally  with  small  step  size  Aq  =  0  01,  until  the  iterations 
stop  to  converge.  In  all  cases,  the  modified  Rayleigh  e(|uation  is  integrated  over  a 
mesh  of  4000  grid  points,  equally  ilislributed  over  the  interval  — l-(^<y<4-(-('i. 
Thus,  the  expected  short  wave  behavior  of  the  additional  instability  modes  is  only 
determined  in  cases  where  their  instability  band  extends  that  of  the  other  modes 

Figures  5-8  show  the  growth  rate  and  phase  speed  of  unstable  modes,  computed 
for  Tr  =  1.  2,4.  and  8.  respectively  For  brevity,  results  obtained  for  the  iiitermediale 
value  S  =  2  are  not  illustrated  In  all  cases,  the  thickness  of  the  heated  layer 
coincides  with  that  of  the  vori'-ily  layer,  tr  =  1.  We  start  with  the  uniform-density 
flow  (Fig  5).  which  IS  later  used  as  reference  to  quantify  the  effects  induced  by  the 
density  variation  Cold  flow  calculations  are  summarized  as  follows: 
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(I)  The  growth  rale  of  the  shear  layer  mode  increases  with  increasing  wake 
deficit.  This  result  is  expected  since  higher  wake  deficits  correspond  to  higher 
vorticity  values. 

The  highest  increase  is  achieved  for  the  smallest  separation  distance,  6  = 
1.  This  IS  not  surprising  and  is  due  to  the  construction  of  the  family  of  velocity- 
profiles  and  the  selection  of  (U',  8)  pairs  so  that  the  same  velocity  extrema  are 
obtained  for  different  separation  distances.  Hence,  smaller  8  values  correspond 
to  higher  vorticity  concentrations. 

(i)  The  wavenumber  of  the  most  unstable  shear  layer  mode  exhibits  a  small 
increase  with  increasing  wake  deficits  at  i  =  1,  but  is  almost  independent 
of  IV  at  higher  separation  distances.  Thus,  the  wake  deficit  does  not  affect 
the  frequency  selection  of  unstable  shear  layer  inode,  which  is  closely  approx¬ 
imated  by  estimates  based  on  the  tank  profile. 

{ij  While  the  growth  rate  of  the  most  amplified  wake  mode  increases  by- 
increasing  W,  Its  wavenumber  of  is  almost  independent  of  H'.  The  growth 
rale  maxima  for  both  the  shear  layer  and  wake  modes  are  reached  for  a  0.6 

(5)  The  phase  speed  of  most  unstable  mode  vanishes  for  the  tank  hyperbolic 
velocity  profile  but  increases  in  the  direction  of  wake  deficit  with  increasing 
VV. 

(6)  The  maximum  growth  rate  of  wake  mode  increases  by  increasing  the  sepa¬ 
ration  between  the  positive  and  negative  vorticity  layers.  However,  in  all  cases 
considered,  the  shear  layer  mode  always  dominates  the  wake  mode.  This  result 
is  best  interpreted  by  focusing  on  the  positive  and  negative  vorticity  layers 
separately,  which  may  be  used  to  distinguish  between  the  symmetric  shear 
layer  and  asymmetric  wake-hke  profiles.  For  the  latter,  the  negative  vorticity 
layer,  whose  inflection  point  is  associated  with  the  unstable  shear  layer  mode, 
has  considerably  higher  strength  and  thus  dominates  wake  component 

The  impact  of  heat  release  on  flowfield  stability  is  examined  in  Figs.  6-8,  which 
show  dispersion  relations  for  Tr  —  2,4,  and  8  The  analysis  is  divided  into  two  sec¬ 
tions;  results  for  the  symmetric  tanft  profile  ate  discussed  first,  and  then  contrasted 
with  corresponding  results  for  the  asymmetric  layer  and  wake  profiles  As  before, 
dispersion  relations  for  the  symmetric  tank  profile  are  lumped  with  tliose  of  the 
asymmetric  shear  layer  inode  with  6=1,  and  identified  by  a  wake  deficit  W  =  0, 
Examination  of  the  dispersion  relations  reveals: 

(1)  By  increasing  T,,  stabilization  of  the  Kelvin-Helinholtz  mode  is  gradu¬ 
ally  achieved.  As  noted  by  McMurtry  et  al.  [9],  who  studied  tlie  stability  of 
a  heated  layer  idealized  by  broken-line  vorticity  and  density  profiles,  heat  re¬ 
lease  inhibits  the  growth  of  the  most  unstable  mode,  whose  amplification  rale 
decreases  to  a  small  fraction  of  the  uniform-density  maximum  as  temperature 
ratio  becomes  high.  The  results  also  indicate  that  the  waveiiuiuliers  of  the 
most  amplified  mode  and  neutrally  stable  modes  decrease  as  Tr  inrrea.ses. 
Thus,  the  density  variation  alters  the  features  of  the  instability  in  such  a  way 
aus  Vo  stabilize  short  wavelength  disturbances  and  collap.se  the  instability  to  a 
thin  baud  affecting  long  wave  perturbations.  On  the  other  hand,  the  vanishing 
phase  speed  property  of  shear  layer  mode  persists. 

(2)  The  stability  properties  of  the  layer  are  significantly  altered  for  Tr  >  4. 
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As  preilicted  in  the  previous  section,  the  computed  results  show  that  the  flow 
field  now  admits  three  unstable  modes.  The  additional  pair  of  unstable  modes 
are  associated  with  the  outer  zeros  of  (rL  'Y  and  are  thus  called  outer  modes 
[8]-  Dispersion  relations  of  these  modes  appear  as  extensions  of  that  of  the 
Kelviu-Helmlioltz  mode. 

(3)  The  outer  modes  have  identical  growth  rates  and  ei(ual  but  opposite 
phcise  speeds.  The  phase  speeds  of  the  neutrally  stable  stalutions  coincide 
with  the  flow  velocity  at  the  outer  inflection  points  Meanwhile,  at  high  Tr 
the  outer  mode  dominates  the  shear  layer  mode  which  is  almost  stabilized 
by  the  density  variation.  The  wavenumber  of  the  most  amplified  outer  modes 
increases  with  increasing  temperature  ratio;  thus,  the  associated  mechanism 
promotes  short -wavelength  instability. 

The  stability  properties  of  asymmetric  shear  layers  and  wake  profiles  differ  sig¬ 
nificantly  from  those  of  the  synunettic  lank  profile,  as  the  heat  release  has  a  less 
significant  impact  than  in  the  former  case  The  stability  results  for  these  profiles  are 
divided  into  two  groups,  according  to  the  separation  distance  between  the  layers  of 
opposite  vorticity  For  a  large  separation  distance  (6  =  3),  the  amplification  curves 
for  the  shear  layer  mode  are  weakly  affected  by  the  presence  of  a  heated  region 
within  the  vorticity  layer.  In  particular,  adl  of  the  stability  properties  of  the  unsta¬ 
ble  modes  are  unaffected,  as  neither  the  instability  bandwidth  nor  the  wave  speed 
of  unstable  modes  are  affected  by  the  temperature  ratio.  This  result  is  expected, 
since  the  density  is  constant  and  equal  to  unity  except  in  the  region  separating  the 
distinct  vorticity  layers.  In  this  region,  the  velocity  profile  is  almost  constant,  so 
that  (p')'  is  closely  approximated  by  the  vorticity  derivative  I'"  in  the  entire  flow. 
Thus,  the  vorticity  maxima  do  not  lie  with  regions  of  high  density  gradient,  and 
results  for  uniform-density  flow  are  recovered. 

For  small  separation  distances,  the  presence  of  a  heated  fluid  layer  results  in 
appreciable  changes  in  the  stability  of  the  flow,  and  has  a  different  influence  on  the 
behavior  of  the  wake  and  shear  layer  modes.  By  increasing  the  temperature  ratio, 
instability  of  the  wake  mode  is  promoted,  as  the  corresponding  amplification  curves 
admit  higher  maxima.  While  the  phase  speeds  of  the  unstable  wake  modes  is  weakly 
affected,  the  corresponding  instability  bandwidth  is  increased,  and  the  wavenumber 
of  the  most  unstable  wtdte  mode  favors  shorter  wavelength  instability.  This  result 
should  be  favorably  contrasted  with  the  results  of  Koochesfahani  and  Frieler  who 
showed  that  the  spatially-developing  wake  mode  may  become  dominant  when  the 
asymmetric  vorticity  layer  is  subjected  to  a  severe  monotonic  density  aiherence  [1]. 
However,  in  all  cases  considered  in  our  study,  the  highest  amplification  rate  of  the 
shear  layer  mode  is  always  considerably  higher  than  for  the  corresponding  wake 
mode.  Therefore,  the  effect  of  heat  release  is  not  expected  to  lead  to  a  qualitative 
change  in  the  behavior  of  the  perturbed  reacting  layer,  which  remains  dominated 
by  the  growth  of  unstable  shear  layer  waves. 

Shear  layer  modes  exhibit  a  different  response  to  the  imposed  density  variation. 
This  response  resembles  the  behavior  of  shear  layer  modes  in  heterogeneous  layers 
for  which  the  density  increases  monotonically  from  one  side  of  the  layer  to  the  other. 
In  both  instances,  density  variation  has  minimal  influence  on  either  the  growth  rate 
of  unstable  Kelvin-Helmholtz  modes,  or  on  their  stability  bandwidth.  However,  it 
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significantly  affects  the  phase  velocity  of  the  waves  In  the  heterogeneous  shear 
layer  [19],  unstable  shew  layer  modes  acquire  an  additional  convective  velocity- 
component  of  the  same  direction  as  the  high  density  stream  The  analogy  between 
the  two  cases  can  be  established  by  inspecting  the  behavior  of  the  density  field  in  the 
neighborhood  of  the  inflection  point  of  the  velocity  profile.  Here,  the  "unstable  shear 
layer  mode"  is  associated  with  the  "upper"  inflection  point  of  the  velocity  profile 
The  density  profile  has  positive  derivative  in  the  neighborhood  of  this  inflection 
point,  so  that  the  density  increases  as  we  move  towards  the  top  free  stream.  Thus, 
we  we  led  to  expect  an  increase  in  the  phase  velocity  of  unstable  shear  layer  mode, 
since  the  top  stream  velocity  is  positive.  This  expectation  is  reflected  in  the  stability 
calculations,  which  show  that  the  phase  velocity  of  the  most  unstable  shew  layer 
mode  increases  with  increasing  temperature  ratios  and  becomes  positive  for  small 
wake  deficits. 

3.2  EFFECT  OF  DENSITY  PROFILE  THICKNESS 

The  stability  results  of  the  previous  section  are  reexamined  for  different  thick¬ 
nesses  of  the  density  profile.  For  laminar  nonpremixed  flames,  the  thickness  of  the 
low-density  zone  depends  on  both  the  thermal  and  mass  diffu.sivities.  Thus,  accurate 
estimates  of  this  thickness  relative  to  the  vorticity  thickness  require  the  solution 
of  the  boundwy  layer  equations  for  reacting  flow,  and  the  results  will  depend  on 
the  Prandtl  and  Lewis  numbers,  and  on  the  details  of  the  chemical  reaction.  In 
general,  the  product  zone  will  be  embedded  within  the  vorticity  layer,  since  near 
unity  Prandtl  and  Lewis  numbers  generally  prevail,  and  initial  conditions  describe 
a  finite  thickness  vorticity  layer  and  a  shwp  interface  separating  the  oxidizer  and 
fuel  streams. 

In  this  study,  however,  such  a  detailed  study  is  replaced  by  the  simplified  ap¬ 
proach  of  considering  different  values  for  the  density  thickness,  which  are  selected 
in  a  wide  parameter  range  in  order  to  cover  most  situations  of  interest  To  this 
end,  the  stability  of  all  the  velocity  profiles  of  Table  I  is  investigated  for  a  variable 
density  field  specified  by  a  fixed  temperature  ratio  Tr  =  4,  which  is  characteristic 
of  a  large  number  of  combustion  appheations.  Meanwhile,  the  density  thickness  is 
gradually  varied;  the  values  cr  =  1.5,0.75,0.5,  and  0.25  are  considered.  The  results 
of  the  computations  are  shown  in  Figs.  9-12,  in  terms  of  the  growth  rate  and  phase 
speed  of  unstable  modts.  As  before,  results  for  the  symmetric  lanh  velocity  profile 
are  lumped  with  those  of  shear  layer  modes  having  ^  =  1. 

The  discussion  of  the  results  distinguishes  between  the  stability  properties  of 
the  shear  layer  and  wake  modes,  and  those  obtained  for  the  symmetric  lanh  profile. 
The  latter  case  is  discussed  first  and  is  summwized  as  follows.  The  behavior  of  the 
tank  shear  layer  mode  is  non-monotonic  with  respect  to  variation  of  the  density 
thickness.  The  results  are  best  interpreted  by  first  considering  the  limiting  case  of 
very  Iwge  density  thickness.  For  this  density  configuration,  the  vorticity  field  lies 
•n  a  zone  of  almost  constant  density,  so  that  the  uniform-density  results  are  recov¬ 
ered.  As  <r  is  decreased  and  becomes  close  to  the  vorticity  thickness,  stabilization 
of  the  shew  layer  mode  is  observed.  The  maximum  growth  rate,  the  most  ampli¬ 
fied  wavenumber  and  the  stability  bandwidth  we  all  reduced.  On  the  other  hand, 
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the  nature  of  the  instability  is  not  altered,  as  all  unstable  modes  have  zero  phase 
velocity.  As  (t  is  further  decreased,  stabilization  of  the  shear  layer  mode  contin¬ 
ues  but  the  flow  field  acquires  an  additional  pair  of  unstable  modes.  Due  to  the 
symmetry  of  the  density  and  velocity  profile,  these  additional  "outer  modes”  have 
identical  instability  bandwidths  and  amplification  curves,  with  equal  magnitude 
but  opposite-sign  phcise  velocities  In  the  range  of  density  thicknesses  considered, 
the  instability  of  the  additional  modes  is  promoted  by  decreasing  a .  as  the  corre¬ 
sponding  most-amplified  modes  admit  higher  growth  rates  and  wavenumbers.  This 
is  accompanied  by  an  increase  in  the  stability  bandwidth  and  the  phase  velocity 
magnitude.  Thus,  tank  shear  layers  which  accommodate  high  heat  release  slow  re¬ 
actions  are  more  susceptible  to  such  short  wavelength  instabilities  However,  we  do 
not  expect  this  mechanism  to  persist  continuously  as  it  is  further  decreased,  since 
in  the  limit  a  [  0  the  deposited  energy  vanishes  so  that  results  for  uniform-density 
flow  should  be  approached. 

The  effect  of  the  density  density  thickness  is  greatly  attenuated  in  the  presence 
of  a  wake  deficit.  In  all  cases,  weak  variations  of  the  stability  properties  of  the  flow 
are  recorded  as  the  value  a  is  altered  The  response  of  the  wake  mode  to  changes  in 
the  density  thickness  cr  is  extremely  weak,  as  the  results  exhibit  almost  insignificant 
changes  in  the  stability  bandwidth,  amplification  curves,  and  phase  relationships 
These  changes  are  of  little  importance  since  the  shear  layer  mode  dominates  the 
initial  evolution  of  the  flow.  The  latter  is  weakly  affected  by  changes  in  the  density 
thickness,  which  result  in  small  modulation  of  the  phase  speed  of  unstable  modes 
As  previously  mentioned,  this  result  is  expected,  based  on  the  similarity  between 
the  behavior  of  the  vorticity  and  density  profiles  around  the  unstable  inflection 
point  and  that  described  in  monotonic  variable-density  shear  layers  [19] 

3.3  VISUALIZATION  OF  UNSTABLE  EIGENFUNCTIONS 

Finally,  the  evolution  of  asymmetric  shear  layers  is  numerically  computed  At¬ 
tention  is  focused  oil  the  late  stages  of  flowfield  evolution,  which  witness  the  for¬ 
mation  of  vortical  structures  due  to  the  non-linear  breaking  of  unstable  waves.  The 
results  are  used  to  examine  the  validity  of  extending  the  trends  established  in  the 
linear  stability  analysis.  This  exercise  is  limited  to  a  visualization  of  the  effect  of 
these  structures  on  the  deformation  of  the  flow  and  the  evolution  of  the  vorticity. 
A  detailed  investigation  of  the  dynamics  of  the  flow  is  not  attempted  in  this  work, 
SIS  detsiiled  reacting  flow  compulations  will  be  discussed  in  a  subsequent  study. 

Numerical  simulation  of  the  variable-density  flowfield  is  performed  using  the 
two-dimensional  transport  element  method.  The  numerical  scheme,  which  belongs 
to  an  adaptive  class  of  Lagrangian  field  methods,  is  based  on  the  discretization 
of  the  vorticity  and  density  gradient  fields  into  a  number  of  transport  elements  of 
finite  overlapping  circular  cores  Accordingly,  the  velocity  field  is  given  by  a  discrete, 
desingularized  convolution  over  the  induced  vorticity  field  of  the  transport  element. 
A  similar  convolution  yields  the  density  field.  Once  the  velocity  field  is  computed  at 
the  element  centers,  a  second-order  predictor  corrector  integration  scheme  is  used 
to  track  their  motion  and  to  advance  the  numerical  solution.  Discrete  vorticity 
and  density  values  evolve  according  to  Eqs.  (4)  and  (2),  respectively.  Meanwhile, 
discrete  density  gradients  are  updated  by  relating  their  evolution  to  the  material 
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deformation  of  the  Lagraiigian  mesh  Details  of  the  formulation  and  construction 
of  this  numerical  scheme,  which  has  been  extensively  employed  in  the  simulation  of 
vtiriable-density  and  reacting  flows,  can  be  found  elsewhere  [13-14,19].  Thus,  only 
a  brief  account  of  the  computational  parameters  used  in  tl.^  calculations  is  given. 

The  visualization  of  unstable  modes  is  performed  using  the  temporal  model  of 
the  vorticity  layers.  In  this  model,  the  velocity  and  density  fields  are  spatially  peri¬ 
odic  in  the  streamwise  direction.  The  periodicity  length.  A,  is  close  to  the  wavelength 
of  the  most  unstable  mode,  as  determined  above.  The  region  of  finite  vorticity  and 
density  gradient  is  initially  discretized  on  a  square  mesh,  having  \ li  and  S S  grid 
points  in  the  cross-stream  and  streamwise  directions,  respectively.  The  core  size  of 
the  transport  elements,  8,  and  the  discrete  values  of  vorticity  and  density  gradient 
are  found  by  minimizing  the  error  between  the  numerical  and  initial  fields  [13]. 

Nine  cases  of  differenitally-heated  layers  and  wakes  are  numerically  investigated. 
The  layers  are  identified  by  values  of  the  temperature  ratio,  the  thickness  of  the 
density  layer,  the  wake  deficit  and  the  displacement  of  the  wake  profile,  which  are 
listed  in  Table  II  alongside  the  wavenumber,  wavelength,  and  growth  rate  of  the 
corresponding  most  unstable  mode.  Table  II  also  shows  the  number  of  computa¬ 
tional  grids  in  the  cross-stream  direction,  NR,  and  the  core  radius  of  the  transport 
elements.  In  all  calculations,  the  time  step  At  =  0  02. 

The  first  three  cases  correspond  to  uniform-density  flow,  and  are  selected  in 
order  to  examine  the  effects  of  the  wake  deficit  and  the  displacement  of  the  wake 
profile.  The  numerical  experiments  are  repeated  by  keeping  the  same  (initial)  vor¬ 
ticity  field  parameters  and  altering  the  initial  density  field  by  letting  Tr  =  4,  and 
6=1.  Thus,  for  these  initial  flow  configurations  (cases  4-6),  the  initial  vorticity 
and  density  thicknesses  coincide.  Finally,  the  dependency  of  the  evolution  of  the 
flow  field  on  the  thickness  of  the  initial  density  profile  is  exanuned  in  cases  7-9. 
These  cases  correspond  to  vorticity  field  configurations  for  which  the  linear  stabil¬ 
ity  analysis  predicts  a  strong  (case  7)  or  weak  (cases  8  and  9)  response  to  changes 
in  6. 

The  computations  are  initialized  by  introducing  a  perturbation  in  the  flow  field 
using  sinewaves  having  the  same  periodicity  wavelength  as  the  computational  do¬ 
main,  A,  and  amplitude  O.OIA.  The  perturbation  is  applied  by  displacing  the  location 
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of  the  iraiisport  elements  m  the  cross-sireeun  direction  according  to  the  sinewaves 
The  calculations  are  extended  in  order  to  observe  the  linear  amplification  of  the 
unstable  modes,  2uid  the  early  non-linear  stages  of  their  evolution  For  cases  1,4,  ^ 

and  7,  this  objective  is  achieved  by  carrying  out  the  computations  until  I  —  12  In 
the  remaining  Ccises,  calculations  are  stopped  at  f  =  9  Results  of  the  simulation  of 
the  flowfield  starting  from  the  flow  coiiflgiirations  of  cases  1-9  are  shown  in  Figs 
13-21,  respectively  The  figures  are  generated  by  plotting  the  location  and  instan¬ 
taneous  position  and  velocity  vector  of  the  transport  elements  The  development  of 
the  vorticity  fi"ld  thus  illustrated  is  discussed  below.  g 

Figure  13  shows  that  the  evolution  uniform-density  shear  layer  initially  de¬ 
scribed  by  the  symmetric  lanh  velocity  profile  does  not  destroy  the  symmetry  of 
the  vorticity  field.  The  vanishing  phase  speed  of  the  instability  wave,  predicted  by 
the  linear  theory  for  the  linear  Kelvin-Hehnholtz  wave,  persists  as  the  waves  un¬ 
dergo  a  non-linear  growth  regime,  and  roll  to  form  a  concentrated  core  of  vorticity 
Detailed  computations  of  similar  flow  fields  have  been  performed  previously  [13], 
and  have  shown  that,  when  pairing  is  disabled,  the  late  stages  are  characterized  by 
a  maturation  of  the  vortex  cores  and  the  continuous  entrainment  of  the  vorticity 
from  the  braids  into  the  cores. 

The  superposition  of  wake  deficit  on  the  synmietric,  monotomc  layer  profile 
results  in  a  significant  departure  from  the  previously  discussed  behavior  The  evo¬ 
lution  of  the  asymmetric  layers  of  Figs.  14  and  15  is  characterized  by  the  finite  wave  § 

speed  of  the  linear  instability  modes,  and  by  the  convective  motion  of  the  vortical 
structures  which  form  following  their  nonlinear  evolution  As  previously  discussed, 
the  direction  of  motion  of  the  linear  instability  waves  can  be  determined  from  the 
the  linear  stability  analysis  and  by  inspection  of  the  initial  velocity  profile.  The 
computations  are  in  agreement  with  the  results  of  the  linear  theory,  which  predicts 
almost  equal  phase  speeds  for  the  most  unstable  modes,  and  a  slightly  larger  growth  ^ 

rate  for  the  layer  with  the  smaller  separation  distance  Alternatively,  this  behavior 
can  be  qualitatively  predicted  by  considering  the  contribution  of  the  "weaker"  vor¬ 
ticity  layer,  whose  induced  flow  field  leads  to  the  motion  of  the  linear  waves  and  of 
the  vortices. 

The  consideration  of  the  effects  of  opposite  regions  of  vorticity  is  easier  for  the 
larger  displacement  parameter  {6  =  3,  Fig.  15).  since  the  corresponding  vortic-  ^ 

ity  profiles  are  formed  of  well-separated  strips  of  opposite  sign  of  vorticity.  The 
computed  results  indicate  that  the  evolution  of  the  flowfield  is  dominated  by  the 
stronger  (negative)  vorticity  layer.  The  latter  appears  to  develop  independently  of 
the  weaker  vorticity  layer,  which  does  not  exhibit  significant  deformation  during  the 
period  of  the  simulation.  This  assessment  no  longer  holds  in  the  case  6  =:  1,  which 
shows  that  both  regions  of  vorticity  deform  simultaneously  (F:g.  14).  This  results  ^ 

in  a  greater  deformation  of  the  vorticity  field  and  the  formation  of  a  substantially 
larger  vortex  core. 

The  impact  of  the  density  variation  on  the  development  of  the  flow  is  illustrated 
in  Figs.  16-18,  which  show  the  evolution  of  shear  layers  characterized  by  the  same 
initial  vorticity  distributions  previously  considered,  and  a  density  profile  having 
Tr  =  4,  and  <7=1.  When  starting  with  the  symmetric  tank  velocity  profile,  the  ^ 

growth  of  Kelvin- Helmholtz  waves  is  significantly  suppressed.  Moreover,  the  non- 


•  •••••• 


354 


linear  wavebreaking  of  these  modes  results  tii  the  formation  of  weak  vortex  cores 
which  are  less  coherent  than  in  the  uniform-flow  case  (Fig  16)  These  results  are  in 
agreement  with  the  predictions  of  the  linear  stability  theory  tuid  with  the  compu¬ 
tations  of  McMurtry  et  al  [9]  who  considered  the  evolution  of  synanetric  reacting 
shear  layers 

Density  variation  effects  are  considerably  less  pronounced  in  the  presence  of  a 
wake  coniponent  For  large  separation  distance,  d  —  i.  comparison  of  the  uniform- 
and  variable-density  results  (Figs.  15  and  18,  respectively)  shows  that  the  presence 
of  a  healed  layer  has  almost  no  influence  on  the  development  of  the  flow  These 
results  extend  those  of  the  linear  stability  theory,  which  predicts  little  changes 
in  the  wavelength,  growth  rate,  and  phaise  speed  of  the  most  unstable  Kelvin- 
Helmholtz  waves  with  the  variation  of  the  density  profile  for  this  initial  votticity 
configuration  On  the  other  hand,  the  interactions  of  the  density  and  vorticity  field 
cannot  be  neglected  for  smadl  separation  distance.  6=1  As  indicated  in  Table  11, 
the  wavelength  of  the  most-unstable  mode  is  decreased  tis  the  temperature  ratio  is 
increased  to  Tr  —  4  .Moreover,  while  the  growth  rate  of  this  mode  is  not  significantly- 
altered  by  the  presence  of  the  heated  layer,  its  phase  speed  is  noticeably  reduced. 
This  observation  also  holds  when  considering  the  cotivectiVe  motion  of  the  vortices 
which  form  following  the  nonlinear  evolution  of  the  unstable  waves  (Figs.  14  and 
17)  However,  in  both  cases,  the  vorticity-density  interactions  do  not  suppress  linear 
growth  and  do  not  inhibit  the  formation  of  large  coherent  vortex  cores 
The  sensitivity  of  the  flow-field  to  the  details  of  the  initial  density  profile  is  examined 
ill  Figs.  19-21.  In  particular,  the  stabilization  of  the  symmetric  layer  by  heal  release 
and  t  he  weak  dependence  of  the  asy  mmet  ric  layer  on  the  presence  of  t  he  heated  layer 
is  investigated.  As  previously  mentioned,  when  the  regions  of  high  vorticity  and 
density  gradient  are  well  separated,  the  initial  development  flow  is  approximated 
by  the  uniform-density  equations,  so  that  layers  characterized  by  a  large  separation 
distance  will  not  be  further  considered.  For  the  symmetric  case.  Fig.  19  indicates 
that  by  decreasing  the  thickness  of  the  density  profile  to  cr  =  0  5,  the  behavior  of 
the  shear  layer  undergoes  an  additional  transition.  The  results  reflect  the  prediction 
of  the  stability  theory  which  indicates  that  in  this  regime  the  most  unstable  mode 
consists  of  a  pair  of  travelling  waves  of  equal  growth  rates,  and  equal  but  opposite 
phase  speeds.  These  waves  amplify  as  they  move  away  from  the  middle  of  the 
computational  domain.  The  rollup  of  the  waves  occurs  as  the  two  trains  of  periodic 
waves  meet  at  the  boundaries  of  the  computational  domain,  and  is  followed  by  the 
formation  of  vortices  whose  cores  are  smaller  and  less  coherent  than  their  uniform- 
density  counterparts.  This  mechanism  differs  significantly  from  the  rollup  of  the 
linear  Kelvin-Hclmlioltz  waves  in  uniform-density  flow,  which,  as  indicated  by  the 
theory,  have  appreciably  larger  growth  rates.  Thus,  the  stabilizing  effects  of  heat 
release  are  expected  to  persist  for  this  profile,  though  the  details  of  the  density- 
distribution  may  lead  to  radical  changes  in  the  development  of  the  flow  and  in  the 
structure  of  the  associated  vorticity  field. 

On  the  other  hand,  the  addition  of  a  wake  deficit  greatly  diminishes  the  "stabilizing" 
impact  of  the  density  variation.  As  predicted  by  the  linear  theory  and  observed 
in  the  computations,  the  early  evolution  of  the  asymmetric  layer  having  a  large 
wake  deficit  is  almost  insensitive  to  the  presence  of  the  heated  layer.  For  such 
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imiial  vorticity  configurations,  the  effects  of  dv-nsity  variation  are  limited  to  a  small 
modulation  of  the  wavelength,  growth  rate,  and  phase  speed  of  the  most  unstable 
mode.  Moreover,  as  illustrated  in  Figs.  14.  17,  20  and  21,  similar  vortical  siructurt's 
are  obtained  aus  a  result  of  the  noii-tinear  wavebreaking  of  the  unstable  modes.  This 
observation  holds  for  all  the  asymmetric  layers  considered,  despite  the  f.ict  that  the 
details  of  the  vorficity  distribution  within  the  vortex  cores  and  in  the  braids  joining 
neighboring  vortices  and  the  convective  motion  of  the  eddies  are  strongly  affected 
by  baroclinic  vorticity  generation  in  the  later  stages  of  evolution  of  the  flow  [13  19] 
4.  Couclusions 

In  this  work,  stability  of  heterogeneous  shear  flows  is  investig.ated  using  linear  stabil¬ 
ity  analysis  and  numerical  simulations  A  large  number  of  initial  flow  ronfiguralions 
which  model  the  development  of  nonpreimxed  reacting  shear  flow  are  analyzed  The 
initial  conditions  are  used  to  exanune  the  effects  of  heat  release  associaieil  with  the 
development  of  a  nonpremixed  flame  on  the  stability  of  asymmetric  shear  layers 
and  wakes.  The  latter  belong  to  a  continuous  family  of  velocity  profiles  which  is 
constructed  by  deforming  the  symmetric  taii/t  shear  layer  profile  using  a  wake  com¬ 
ponent  characterized  by  characterized  by  a  spread.  and  a  velocity  deficit.  VT 
Density  variation  is  used  to  model  the  effects  of  heat  release  1  he  initial  density 
distribution  corresponds  to  a  temperature  spike  and  is  described  by  a  thickness,  <t, 
and  a  temperature  ratio,  Tr  Stability  curves  are  obtained  m  this  four-dimensional 
parameter  space  and  numerically  visualized  using  the  transport  element  method. 
The  numerical  simulations  are  extended  into  the  nonlinear  stages  of  flowfield  evo¬ 
lution  in  order  to  examine  the  validity  of  extrapolating  the  linear  stability  results 
Stability  of  the  lanh  reacting  shear  layer  exhibits  strong  sensitivity  to  the  details 
of  the  density  distribution.  When  the  density  and  vorticity  thicknesses  are  close, 
stabilization  of  the  shear  layer  occurs  as  the  temperature  ratio  increases  This  effect 
is  manifested  by  a  sharp  decrease  in  the  instability  growth  rate,  the  nonlinear  evo¬ 
lution  of  unstable  eigenfunctions  yields  weaker  less-coherent  vortex  structures  than 
those  observed  in  uniform-density  flow  By  decreasing  the  thickness  of  the  density 
profile,  additional  inflection  points  and  instability  modes  are  observed  These  insta¬ 
bility  modes  dominate  the  shear  layer  mode  which  is  almost  stabilized  by  the  heat 
release.  However,  while  the  associated  growth  rates  increase  w  ith  decreasing  density 
thickness  and  appear  to  approach  the  maximum  growth- rate  values  obtained  for 
uniform-density  flow,  the  nonlinear  evolution  of  these  instability  modes  does  not 
result  in  the  formation  of  concentrated  vortices  or  in  substantial  deformation  of  the 
flow  Thus,  for  the  taii/i  profile,  heat  release  tends  to  stalnlize  of  the  flow 

Addition  of  a  wake  deficit  significantly  alters  the  stability  of  the  flow.  In  the 
parameter  range  considered,  two  flow  configurations  are  distiiigiiished  When  the 
zones  of  high  density  gradients  and  vorticity  magnitude  are  well  separated,  heat 
release  and  density  variation  have  almost  no  impact  on  flowfield  stability.  In  this 
case,  the  growth  rate  and  phase  speed  of  unstable  waves  are  weakly  affected,  and 
their  non-linear  evolution  results  in  the  formation  of  large  concentrated  asymmetric 
vortices.  On  the  other  hand,  when  regions  of  high  vorticity  and  density  gradients 
are  close  or  coincide,  heat  release  has  a  weak  effect  on  the  development  of  the  flow 
While  the  growth  rates  of  unstable  waves  are  almost  unaffected,  their  phase  speeds 
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depend  on  the  details  of  the  density  variation  This  mechanism  is  also  observed  in 
the  later  stages  of  evolution  of  the  flow,  which  indicate  that  the  convective  motion 
of  the  large  vortices  is  modulated  in  the  same  manner  as  the  phase  velocity  of  the  ft 

corresponding  linear  waves  Thus,  heal  release  is  not  expected  to  inhibit  the  growth 
of  unstable  modes  in  asymmetric  lajers  and  wakes,  but  may  however  influence  the 
global  features  of  the  flow 

The  correspondence  between  computed  stability  results  and  initial  flow  config¬ 
urations  indicates  that  the  effects  of  heat  release  may  be  intuitively  predicted  by 
simple  examination  of  the  behavior  of  the  density  profile  in  the  neighborhood  of  the  ft 
inflection  point  of  the  velocity  profile  If  the  density  profile  reaches  a  mniitnum  in 
the  neighborhood  of  the  inflection  point  of  the  velocity  profile,  stabilization  of  the 
corresponding  >>!>stable  mode  is  expected  If  the  density  gradient  does  not  vanish  in 
this  neighborhood,  we  expect  a  minor  variation  in  the  growth  rate  of  the  unstable 
mode,  and  a  modulation  of  its  phase  speed  which  depends  on  the  sign  of  the  density 
gradient  In  this  case,  the  results  exhibit  analogous  trends  to  those  established  for  ft 
symmetric  shear  layers  separating  streams  of  unequal  density,  where  unstable  waves 
are  observed  to  acquire  a  streaniwise  convection  component  in  the  direction  of  the 
high-density  stream  (13]. 

The  large  predicted  differences  in  the  stability  properties  for  syiiuiietric  and 
asymmetric  reacting  layers  lead  us  to  expect  significant  sensitivity  of  developing 
reacting  shear  layers  to  initial  disturbances  In  particular,  if  unstable  modes  are  ft 
excited  at  short  distances  downstream  of  the  splitter  plate,  i  e  before  viscous  dif¬ 
fusion  leads  to  the  destruction  of  the  wake  deficit  associated  with  the  merger  of  the 
Blasius  profiles,  significant  stabilization  by  heat  release  is  not  expected  Otherwise, 
a  sharp  decrease  in  shear  layer  growth,  mixing  and  burning  efficiency  is  anticipated 

We  finally  note  that,  in  the  parameter  range  considered,  shear  layer  modes  were 
always  found  to  dominate  the  wake  modes.  Thus,  heat  release  and  density  variation  ft 
are  not  expected  to  produce  a  significant  change  in  the  shape  of  instability  [»).  Our 
results  should  be  contrasted  with  those  of  Koochesfahani  and  Frieler  [1]  who  sitowed 
that  for  large  wake  deficits  and  density  differences,  the  initial  development  of  the 
layer  can  be  dominated  by  the  amplification  of  unstable  wake  modes  A  stability 
analysis  in  this  flow  regime  was  omitted,  since  the  corresponding  flow  configurations 
are  not  representative  of  typical  reacting  shear  flow  applications.  Extension  of  the  ( 
parameter  range  to  include  these  initial  conditions  is  postponed  to  a  subsequent 
study  which  will  focus  on  direct  simulation  of  the  developing  reacting  flow 
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Figure  1 .  Schematic  sketch  of  a  shear  layer  ( left )  ami  a  blutf-Lody  wake  flow  ( right ) 
The  self-similar  monotomc  sliear  layer  profile  is  shown  in  curve  (a),  while  curves 
(h)  ami  (c)  respectively  illustrate  asyiiimetric  shear  layer  ami  wake  profiles  The 
top  (bottom)  row  corresponds  to  ^  =  1  (#  =  3) 


Figure  2.  Velocity  profiles  for  the  family  of  shear  layers  givi  n  by  F,(|  (lU).  with 
(H  ,  d)  combinations  of  Table  I  The  plots  are  respectively  arranged  from  top  for 
increasing  wake  increasing  thickness,  f  =  1,2.  and  3  The  tniih  profile  is  shown  with 
the  =  1  subcollection.  The  top  (bottom)  row  corresponds  to  ('>  =  1  (fi  =  3) 
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Figure  3.  Profiles  of  f"  +  for  layers  described  by  <t  =  1 ,  ^  =  1  and  (a)  Tr  =  1 

(uniform-density  flow);  (b)  Tr  =  2;  (c)  Tr  =  4;  and  (d)  Tr  =  8.  Similar  behavior  is 
obtained  for  the  larger  separation  disctances,  6  =  2  and  3.  The  top  (bottom)  row 
corresponds  to  ^  =  1  (6  =  3). 


Figure  4.  Behavior  of  U"  +  p'U'/p  for  the  velocity  profiles  of  Fig.  2,  with  density 
distribution  given  by:  Tr  =  4,  and  (a)  a  =  1.5;  (b)  cr  =  0.75;  (c)  <t  =  0.5;  and  (d) 
ff  =  0.25.  The  top  (bottom)  row  corresponds  to  i  =  1  (6  =  3). 


Figure  5a.  Growth  rate  vs.  wavenumber  of  the  shear  layer  (left)  and  wake  (right) 
modes  for  the  velocity  profiles  of  Fig.  2  and  uniform-density  profile  T,  =  1.  The 
top  (bottom)  row  corresponds  to  6  =  1  (i  =  3). 


Figure  5b.  Phase  speed  vs.  wavenumber  for  the  shear  layer  (left)  and  wake  (right) 
modes  for  the  velocity  profiles  of  Fig.  2  and  uniform-density  profile  =  1.  The 
top  (bottom)  row  corresponds  to  ^  =  1  (t  =  3). 
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Figure  6a.  Growth  rate  vs.  wavenumber  of  the  shear  layer  (left)  and  wake  (right) 
modes  for  the  velocity  profiles  of  Fig.  2  and  density  profile  given  by  T,  =  2  and 
<T  =  1,  The  top  (bottom)  row  corresponds  to  6  =  1  (6  =  3). 


Figure  6b.  Phase  speed  vs.  wavenumber  of  the  shear  layer  (left)  and  wake  (right) 
modes  for  the  velocity  profiles  of  Fig,  2  and  density  profile  given  by  Tr  =  2  and 
(T  =  1.  The  top  (bottom)  row  corresponds  to  6  =  1  (6  =  3). 


Figure  7b.  Phase  speed  vs,  wavenumber  of  the  shear  layer  (left)  and  wake  (right) 
modes  for  the  velocity  profiles  of  Fig.  2  and  density  profile  given  by  Tr  =  4  and 
cr  =  1.  The  top  (bottom)  row  corresponds  to  6  =  1  (^  =  3). 
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Figure  9a.  Growth  rate  vs.  wavenumber  of  the  shear  layer  (left)  and  wake  (right) 
modes  for  the  velocity  profiles  of  Fig.  2  and  density  profile  given  by  =  4  and 
<r  =  1.5.  The  top  (bottom)  row  corresponds  to  i  =  1  (iS  =  3). 


Figure  9b.  Phase  speed  vs.  wavenumber  of  the  shear  layer  (left)  and  wake  (right) 
modes  for  the  velocity  profiles  of  Fig.  2  and  density  profile  given  by  Tr  =  4  and 
(T  =  1.5.  The  top  (bottom)  row  corresponds  to  i  =  1  (6  =  3). 


Figure  10b.  Phase  speed  vs.  wavenumber  of  the  shear  layer  (left)  and  wake  (right) 
modes  for  the  velocity  profiles  of  Fig.  2  and  density  profile  given  by  Tr  =  4  and 
a  =  0.75. 
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Figure  13.  Evolution  of  the  uniform-density  shear  layer  with  velocity  profile  given  by 
W  =  i  =  0,  illustrated  in  terms  of  the  vortex  elements.  The  plots  are  generated  by 
drawing  the  location  and  instantaneous  velocity  vector  of  the  elements  at  t  =  4.8, 
and  12. 
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Figure  14.  Vortex  element  representation  of  the  uniform-density  shear  layer  with  §  ^ 

velocity  profile  given  by  IV  =  1.2  and  6  =  1,  at  t  =  3,6,  and  9  The  plots  are 
generated  as  in  Fig.  13. 
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Figure  15.  Vortex  eiemertt  representation  of  tlie  uniform-density  shear  layer  with 
velocity  profile  given  by  W'  =  0  82149  and  #  =  3,  at  t  =  3,6.  and  9.  The  plots  are 
generated  as  in  Fig.  13. 
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Figure  16.  Vortex  element  ^“presentation  of  the  shear  layer  with  velocity  profile 
given  by  14'  =  ^  =  0  and  density  profile  having  TV  =  4  and  =  1.  The  plots  are 
generated  as  in  Fig.  13. 


Figure  17  Vortex  element  representation  of  the  shear  layer  with  velocity  profile 
given  by  U’  =  1.2  and  b  =  1  and  density  profile  having  Tr  —  -1  and  ir  =  1  The 
plots  are  generated  as  in  Fig.  13. 


Figure  18.  Vortex  element  representation  of  the  shear  layer  with  velocity  profile 
given  by  =  0.8214')  and  b  =  Z  and  density  profile  having  Tr  -  A  and  a  =  1.  The 
plots  are  generated  as  in  Fig.  13. 


Figure  19.  Vortex  element  representation  of  the  shear  layer  with  velocity  profile 
given  by  IV  =  6  =  0  auid  density  profile  having  Tr  =  4  and  a  —  0  5  The  plots  are 
generated  as  in  Fig.  13. 


Figure  20.  Vortex  element  representation  of  the  shear  layer  with  velocity  profile 
given  by  W  =  12  and  ^  =  1  and  density  profile  having  Tr  =  4  and  cr  =  0,5.  The 
plots  are  generated  as  in  Fig  13. 
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Figure  21.  Vortex  element  representation  of  the  shear  layer  with  velocity  profile 
given  by  IV  =  1.2  and  6  =  1  and  density  profile  having  Ti-  =  4  and  a  =  15,  The 
plots  ate  generated  as  in  Fig.  13. 
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